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Abstract 


The desire to create more complex visual scenes in modern flight simulators 
outpaces recent increases in processor speed. As a result, simulation transport delay 
remains a problem. New approaches for compensating the transport delay in a flight 
simulator have been developed and are presented in this report. The lead/lag filter, the 
McFarland compensator and the Sobiski/Cardullo state space filter are three prominent 
compensators. The lead/lag filter provides some phase lead, while introducing significant 
gain distortion in the same frequency interval. The McFarland predictor can compensate 
for much longer delay and cause smaller gain error in low frequencies than the lead/lag 
filter, but the gain distortion beyond the design frequency interval is still significant, and 
it also causes large spikes in prediction. Though, theoretically, the Sobiski/Cardullo 
predictor, a state space filter, can compensate the longest delay with the least gain 
distortion among the three, it has remained in laboratory use due to several limitations. 

The first novel compensator is an adaptive predictor that makes use of the Kalman 
filter algorithm in a unique manner. In this manner the predictor can accurately provide 
the desired amount of prediction, while significantly reducing the large spikes caused by 
the McFarland predictor. Among several simplified online adaptive predictors, this report 
illustrates mathematically why the stochastic approximation algorithm achieves the best 
compensation results. A second novel approach employed a reference aircraft dynamics 
model to implement a state space predictor on a flight simulator. The practical 
implementation formed the filter state vector from the operator’s control input and the 
aircraft states. The relationship between the reference model and the compensator 
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performance was investigated in great detail, and the best performing reference model 
was selected for implementation in the final tests. 

Theoretical analyses of data from offline simulations with time delay 
compensation show that both novel predictors effectively suppress the large spikes 
caused by the McFarland compensator. The phase errors of the three predictors are not 
significant. The adaptive predictor yields greater gain errors than the McFarland predictor 
for short delays (96 and 138 ms), but shows smaller errors for long delays (186 and 282 
ms). The advantage of the adaptive predictor becomes more obvious for a longer time 
delay. Conversely, the state space predictor results in substantially smaller gain error than 
the other two predictors for all four delay cases. 
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Preface 


This report is the first of two NASA contractor reports documenting the research 
on flight simulator transport delay compensation, undertaken in the Man-machine 
Systems Research Laboratory at the State University of New York at Binghamton and 
supported by the NASA Langley Research Center, in Hampton, Virginia. Loosely 
speaking, the two reports cover the theoretical research and the experimental testing of 
the research, respectively. 

This report begins with a theoretical investigation of the effects of pure time delay 
on a control system consisting of an aerodynamic model, a pilot model and the Pade 
approximation of time delay. It then summarizes the literature study of transport delay 
causes in, and effects on, a flight simulator. This report continues with the introduction of 
three existing transport delay compensators — the lead/lag filter, the McFarland predictor 
and the Sobiski/Cardullo predictor, including intensive analyses of the strengths and 
limitations of each compensator. After a brief description of an expedient algorithm, 
designed to reduce the large spikes by the McFarland predictor, it presents the main body 
of research, i.e., development of two novel compensators. This report then thoroughly 
develops the adaptive predictor and the state space predictor. The adaptive predictor is a 
special Kalman filter that recursively updates the coefficients so that accurate prediction 
can be achieved. Among several versions of the adaptive algorithms, the Stochastic 
Approximation algorithm is mathematically demonstrated to achieve the best 
compensation results. The state space predictor makes use of the state transition matrix 
and its integral of a reference aircraft model. Several aircraft models were tested and the 
landing model of a large commercial transport in pitch achieved the best compensation 
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results as a reference model. By simplifying the state space predictor, the relationship 
between the compensation quality and the reference model was intensively investigated. 
Offline compensation results are presented to compare the McFarland predictor and the 
two novel predictors. The final part of the first report draws conclusions and suggests 
possible future research. 

The second in the series, i.e., NASA CR 2007-21 5096 1 is presented in three parts: 
transport delay measurement in the NASA Langley Research Center’s Visual Motion 
Simulator (VMS), piloted testing of the time delay compensators, and conclusions. The 
time delay measurement was conducted to verily the actual transport delay prior to the 
application of compensation in the final piloted tests. The average transport delay from 
the pilot control input to the visual display update was measured to be 90 ms. The second 
part of the report treats the final piloted experiment design, added time delay, test 
subjects, compensators, data collection, and evaluation metrics. It then presents the 
results of the final piloted tests in terms of performance errors, task load index, handling 
quality and power spectral density of the pilot controls. The final part of the report draws 
conclusions on the delay measurement and piloted simulation tests, and includes 
suggestions for future research. The appendices of the report include resultant graphs of 
all 13 pilots in terms of the four metrics, and the source code and flowcharts of some of 
the algorithms used in the research. 
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1. Introduction 


1.1. Transport Delay in Vehicle Simulation 

The transport delay in a vehicle simulator is the time elapsed from an operator’s 
control input until an appropriate stimulus is presented to the operator by the associated 
hardware - . In a real vehicle, the transport delay is negligible because the vehicle responds 
to the operator command almost instaneously. Unfortunately, this is not true for a vehicle 
simulator. As an example, Fig. 1.1 shows a motorcycle simulator. Unlike a driver on a 
real moving motorcycle who directly feels the motion of the motorcycle relative to the 
street, the driver on this simulator perceives the motion primarily based on the visual 
display showing the movement of the road and the surroundings. The time it takes for the 
simulator computers to generate a new visual image on the screen based on the operator’s 
control input is the transport delay. To illustrate the sources of the transport delay, Fig. 1.2 
shows the architecture of an ordinary vehicle simulator with a visual system. 

The transport delay comes primarily from three sources: sampling delay, 
processing time and data transfer time. Sampling delay results because the simulator 
dynamics computer only samples the operator’s control input at the beginning of each 
computation frame whereas the actual control input arrives stochastically. Therefore the 
change of input between two consecutive sampling events is delayed. It may be as long as 
almost a full frame, or as short as zero, but the average of the sampling uncertainty is a 
half frame. The processing time consists of two parts — the time taken by the dynamics 
computer to calculate the vehicle states from the sampled operator’s control input, and 
the time for the computers in the visual system to prepare the visual image. The 



processing time usually dominates the total transport delay. Data transfer time is the time 
it takes for the visual system to receive the updated vehicle state computed by the 
dynamics computer. If the update rates of the vehicle dynamics computer and the visual 
system are not equal and the latter is not an integer multiple of the former, 
communication asynchrony occurs which results in additional delay. If the transfers are 
asynchronous, the data transfer delay affects the sampling delay. As long as the transfer 
time is less than the sampling interval (i.e., the frame length), transfer time may be 
considered the same as processing time. Although the simulator time delay consists of 
several components from different subsystems, the origin makes no difference to the 
operator, who only feels the total effect. 



Fig. 1.1. A motorcycle simulator 

In Fig. 1.2, the sampling delay occurs between the hand and the plant, the 
processing delay occurs in the plant and between the output and the display, and the delay 
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due to data transfer may arise if there is a difference in update rate between the plant and 
the display system. 



Fig. 1.2. Architecture of a vehicle simulator with a visual system 

If the overall delay reaches a noticeable level, when the operator tries to perfonn a 
task, say a left turn, she will see insufficient response from the display relative to her 
expectations; hence the operator’s cognitive control logic causes her to maneuver further 
until the expected display is observed; but because of the delay, the display will show the 
operator that she has already over controlled, resulting in a compensation or a 
modification, and so on. The resulting locus of the motorcycle positions might resemble 
the dashed curve in Fig. 1.3. 

This example demonstrates that one of the immediate effects of long transport 
delay is Pilot Induced Oscillation (PIO). As the time delay gets longer, the oscillation is 
expected to be more severe — with a larger magnitude and a slower decay, which may 
even become unstable. In other words, time delay makes the system’s response slower 
and undermines the system stability. As a result, the virtual vehicle is harder to control 
with time delay, indicating that the operator’s perception of the handling quality becomes 
worse, and the control workload is increased. Using Fig 1.3, it is easy to visualize the 
degradation in performance by comparing the actual trajectory, represented by the dashed 
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curve with the ideal trajectory, represented by the solid curve. In summary, the following 
problems caused by the time delay are expected: 

1) The man-machine system performance is degraded; 

2) The operator’s control workload is increased due to over control and modification; 

3) The operator’s assessment of the handling quality of the system is diminished. 



Fig. 1.3. Ideal turn of a real motorcycle versus an actual turn in a simulator 

In the frequency domain, time delay shifts the time-line of the simulated vehicle 
to the right, with respect to the response of the real vehicle, causing a phase lag in the 
simulation system. This phase lag is proportional to the frequency components of the 
operator’s control input. The phase lag at the system crossover frequency decreases the 
system phase margin; it also contributes to the PIO and undermines system stability. To 
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restore the system phase margin, the operator tries to increase the control gain or lower 
the system crossover frequency, resulting in an increase in control workload, and 
degrading the handling quality. The frequency analysis agrees with the time domain 
analysis. 

The literature supports the above analyses of the effects of transport delay on a 
man-in-the-loop flight simulator system. Several metrics indicate that transport delay 
degrades the man-machine system performance. Transport delay increases the system 
Root Mean Square Error (RMSE) associated with various tasks (Riccio, et, al, Bailey, et 
al); the Power Spectral Density (PSD) analyses of the operator controls demonstrate that 
the time delay makes the operator’s workload increase, especially in the high frequencies 
(Middendorf, et al, Guo, et al); the Cooper-Harper Rating (CHR) also shows that the 
operator’s handling quality assessment is affected by the delay (Cooper and Harris). 
Large transport delays may also induce simulator sickness (Zaychik, et al). (The literature 
study of the time delay effects is elaborated in Chapter 2.) 

1.2. Delay Compensation 

Because the impact is undesirable in flight simulations, simulator transport delay 
must be minimized in order to reduce its effects. If, after minimization, the transport 
delay still exceeds the tolerable threshold for maintaining desirable simulator 
performance, algorithms to compensate for the delay should be employed. Delay 
compensation usually makes use of prediction of the aircraft states before they are output 
to the cueing channels. This is illustrated in Fig. 1 .4, where, in the small plot to the right 
of the predictor block, the black dashed curve is the predicted aircraft state. Images based 
on the predicted aircraft state can be used to offset the transport delay in the visual system. 
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The purpose of prediction is to restore the phase margin, which would be reduced by the 
transport delay. 

Prediction is achieved by making use of past and current system information, 
including the aircraft displacement, velocity and acceleration, the operator control input, 
the dynamic model, and so on. What information is used for the prediction and how to 
use it lead to various ways of designing the compensator. The most prominent three 
compensators are the lead/lag filter, the McFarland predictor and the Sobiski/Cardullo 
predictor. 



Fig. 1.4. Delay compensation based on prediction 

The lead/lag filter had long been used in industry before Ricard/Harris introduced 
it to the flight simulator to compensate for the transport delay. Having a single pole and a 
single zero, the lead/lag filter provides some phase lead in a certain frequency range 
while introducing gain distortion. In order to properly design a lead/lag filter, the 
designer must determine the pole, the zero and the gain appropriate for the transport 
delay to be compensated such that the phase lead and gain distortion are well balanced. 
Both Ricard/Harris and Crane proposed methods of designing the lead/lag filter, and they 
tested the compensation against the performance of piloted simulations. Though both 
methods show some advantages, the lead/lag compensator has been replaced by other 
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more powerful predictive filters, primarily due to its limited ability to provide phase lead 
and the undesirable introduction of significant gain distortion. McFarland developed a 
discrete filter, which extrapolates the future aircraft displacement, from three consecutive 
iterations of velocity. This special integration algorithm is a type of finite impulse 
response filter because it only has poles at the origin. Because the current prediction does 
not involve the past predictions, the prediction error is not passed to the next iteration, 
and therefore, there is no error accumulation. The large gain distortion, which would be 
present when using the lead/lag filter is significantly reduced while phase lead is 
substantially increased. The challenge in designing this type of filter is to detennine the 
three coefficients that multiply the three steps of velocity. McFarland introduced a 
method known as sinusoidal tuning, which makes use of boundary conditions of the so- 
called “pass band”. The pass band is defined to be the primary frequency band for most 
pilot operations. While the McFarland filter works well within this pass band, the gain 
distortion and phase lead deficiency are significant, and the gain distortion leads to very 
disturbing spikes in the prediction. The spikes originate from the constant coefficients, 
which were determined using the sinusoidal tuning, and are not adjusted during the 
simulation. 

The Sobiski/Cardullo predictor is the first state space filter used for compensating 
the transport delay in a flight simulation. It was derived from the solution of a linear 
time-invariant (LTI) differential equation in state space fonnat. By using more 
information, in each iteration of prediction, theoretically the full order Sobiski/Cardullo 
filter should achieve better compensation (sufficient phase lead and less gain error) than 
prior techniques, provided that the aerodynamics are also LTI and known. 
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However, there are three practical problems which prevent the widespread 
application of the Sobiski/Cardullo predictor to flight simulators. First, most modern 
flight simulators include complex, nonlinear, time variant aircraft models. Second, when 
using the Sobiski/Cardullo predictor, the extrapolated state is only valid if the operator’s 
control input is piece- wise constant, sinusoidal, or exponentially decaying, etc. Third, the 
matrix operations used to implement the Sobiski/Cardullo predictor make it 
computationally intensive, and simplifying the algorithm would make it more practical. 

1.3. Scope of Research 

This is a comprehensive study of the transport delay in a vehicle simulator, from 
its sources, to its effects, measurement and compensation. In Chapter 2, a theoretical 
analysis of a pure time delay — its effects on a control system in both the time and 
frequency domain — is presented. The second part of Chapter 2 is a summary of a 
literature study on the causes and effects of the transport delay in a flight simulator. 

Chapter 3 describes the three prominent compensation techniques, the lead/lag, 
McFarland and Sobiski/Cardullo filters, which were briefly introduced in this chapter 
(Section 1.2), in much more detail. The basic principles, the formulation, and the 
advantages and disadvantages of each filter will be presented in this chapter. Analyses of 
these filters in both the time and frequency domains are also presented. 

The equation proposed by Crane for positioning the pole of the lead/lag filter has 
been revised, and a filter designed with the revised equation shows obvious improvement 
over those designed using Crane’s original equation. The revision is introduced in 
Chapter 3. 



Chapter 4 introduces two novel predictors for compensating transport delay. 
Section 4.1 presents a simple spike reduction algorithm to alleviate the gain distortion 
caused by the McFarland compensator. In section 4.2, a new adaptive predictor is 
introduced. The Kalman estimator, which is an online recursive least squares method has 
been adopted to design the coefficients of a predictive compensator, which also uses three 
consecutive velocities to extrapolate, similar to the McFarland compensator. While 
simplifying the Kalman filter algorithm, a forgetting factor, the Kaczmarz algorithm, the 
stochastic approximation algorithm and the Least Mean Squares algorithm are introduced. 
This section also mathematically demonstrates why the stochastic approximation 
algorithm stands out above all the other adaptive algorithms in compensating the 
transport delay. 

A second novel approach employes a reference aircraft dynamic model to 
implement a state space predictor for use on a flight simulator. The practical 
implementation formed the filter state vector from the operator’s control input and the 
aircraft states. Among several reference models tested, the landing model of a large 
commercial transport in the pitch axis achieves the best compensation result, and was 
selected for the final piloted tests. The relationship between the reference model and the 
compensation performance is also investigated in detail in chapter 4. 

Theoretical analysis of the two novel compensators is the main topic of Chapter 5. 
It covers an evaluation and comparison of the errors caused by a compensator with 
different predictors, and includes a sensitivity analysis, which demonstrates how the 
compensation errors change as the time delay increases. The chapter begins by defining 
two error metrics, and then compares the compensation errors in terms of the two error 
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metrics using offline tests (see Chapter 5). The chapter concludes with a sensitivity 
analysis. 
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2. Background Information 


2.1. Theoretical Description of Time Delay 

A pure time delay simply causes a signal to shift right on the time line. As an 
example, Fig. 2.1 shows a sinusoid signal ( y(t) = sin(cot )) and the resulting signal when it 

is delayed by t d . The delayed signal, y d (t) = sin{cot-cot d ) , has a phase lag of cot d with 

respect to the initial signal. This example demonstrates that time delay can be described 
in both time domain and frequency domain. The transfonn of a time delay t d between 
the time domain and the frequency domain is given in Eq. (2.1). 

i 
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Fig. 2.1. A sinusoid signal and its delayed result by t d 

f(t-t d )^e~ ja,J F{cot d ) (2.1) 

If the signal is continuous, the delayed signal is simply given by 


Time Delay and Phase Lag 
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yd( t ) = y( t ~ t i) 


( 2 . 2 ) 


And for a continuous signal, the relationship between it and its delayed partner in the 
frequency domain is given by the Laplace transfer function 

M*) 


Y(s) 


= e 


(2.3) 


Because e is nonlinear, it is usually approximated by the 2 nd -order Pade approximation 


2 6 12 

s' s + 

Ki 

2 6 12 

S' H 5 + ^r 


(2.4) 


In a discrete system, if the time delay is an integer multiple of the frame time T , 
say t d = n d T , where n d is an integer, the counterparts of Eq. (2.2) and (2.3) are simply 


given by Eq. (2.5) and (2.6) respectively 

y d (k) = y{k-n d ) 

r(z) 


(2.5) 


(2.6) 


However, if the ratio r = ^ is not an integer, these relationships become much 


more complicated. Substituting the trapezoidal integration, ~ = —^~ into Eq. (2.4), the 

s 2 z-l 

Pade approximation becomes 


Y d (z) fr+ftz-'+z- 2 
7(z) 1 + a x z~ x + oc 2 z~ 2 


(2.7) 


where 


oc 2 /? 0 


3r 2 — 3r + 1 
3r 2 +3r + 1 


and a x = /) 


6 r 2 - 2 
3r 2 + 3r + 1 


And the difference equation 


corresponding to Eq. (2.5) is given by 
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y d i k ) = -a, y d (k - 1) - a 2 y d (k - 2) + 0 a y (k) + y (k - 1) + 0 2 y (k - 2) 


(2.8) 


It follows from Eq. (2.3) that the transfer function of a pure time delay has unity 
magnitude at all frequencies, but has negative phase angle, calculated by 

</>,,= t d co (2.9) 

This can be verified by using the Bode diagram of a time delay as shown in Fig. 2.2, in 
which both the exact calculation and the 2 nd -order Pade approximation are plotted. When 
a time delay is added to an open loop system, it only delays the output without causing 
any gain distortion. Take Fig. 2.1, for example, if a time delay of t d is applied to the solid 
sinusoid signal, it moves the curve to the right and becomes the dashed curve. 
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Bode Diagram of a Pure Delay 
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Fig. 2.2. Bode diagram of a pure delay: exact calculation and Pade approximation 

However, if time delay is introduced in a closed loop system, the system output is 
shifted to the right, and the gain changes, because the system feedback is also delayed. 
Delayed feedback makes the system sluggish so that it becomes more oscillatory and its 
stability is undermined. If the delay is sufficiently large, the system could become 
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unstable. This can be illustrated by modeling a flight simulation task, as shown in Fig. 2.3. 
The pilot model, given by Eq. (2.10), matches a lateral control task performed with a rate 
controller cascading a delay term representing the neuromuscular and cognitive time 
delay, which were lumped into the predictor. The aircraft model, given by Eq. (2.11) 
represents the change in the roll angle per unit of deflection of the control stick, at a flight 
condition of 430 knots airspeed and 30,000 feet altitude. The time delay block refers to 
the artificially inserted transport delay (denoted by t d ) represented by the 2 nd -order Pade 
approximation (Eq. (2.12)). Three values of the artificial delay were tested: 0, 200, or 400 
ms, and the closed-loop step responses and the open loop frequency responses of these 
three cases are given in Fig. 2.4 and Fig. 2.5. 



Fig. 2.3. Block diagram of a simulation with a man-in-the-loop control 
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Bode Diagrams of a Control System with Different Delays 
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Fig. 2.4. Bode Diagrams of a closed loop system with different delays 
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Step Responses of a Control System with Different Delays 



Fig. 2.5. Step responses of a closed loop system with different delays 

From the Bode diagrams of the open loop system, it is appearent that the 
magnitude is not changed by the time delay, but the phase angle is decreased, which 
agrees with the unity magnitude and negative phase angle properties of the delay. 
According to Eq. (2.9), the phase margin decrease at the crossover frequency is 
proportional to the amount of time delay 

</>PM=t d a c ( 2 . 12 ) 

For a 200 ms delay, the phase margin is reduced considerably but is still positive, 
which means the system is still stable, yet becomes more oscillatory. When the delay is 
400 ms, the phase margin is negative, which indicates an unstable system. The step 
responses of the closed-loop system confirm the frequency domain analysis (Fig. 2.5). 
Note that although for 200 ms delay the step response has more oscillations, it still 
converges to the same steady state value. 
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The impact of time delay on the man-machine system is not only dependent upon 
the amount of the delay, but also the system dynamics. Specifically, the system 
bandwidth or crossover frequency impacts the effects of time delay on the system. To 
illustrate, a second system is created by varying the real part of the dominant poles of the 
aircraft model given in Eq. (2.11) (i.e., by changing the natural frequency from 1.8788 
rad/s to 1.95 rad/s) with the operator model unchanged. The bandwidth of the new 
system (System II) is larger than that of the original system (System I). For comparison, 
some properties of the systems are listed in Table 2.1. The table shows that time delay 
decreases the system closed-loop bandwidth, and the longer the time delay, the larger the 
decrease in close-loop bandwidth. In addition, the system with higher bandwidth tends to 
suffer faster closed-loop bandwidth reduction with time delay. 

Table 2.1. Properties of two dynamic systems 


Properties 

System I 

System II 

Damping ratio 

0.2376 

0.2376 

Natural frequency (rad/s) 

1.8788 

1.9500 

Crossover frequency (rad/s) 

2.19 

2.54 

Closed-loop bandwidth (rad/s) 

4.3238 

4.4864 

Closed-loop bandwidth (rad/s) 
with 150 ms delay 

3.9131 

4.0338 

Closed-loop bandwidth (rad/s) 
with 300 ms delay 

3.4760 

3.5801 

Phase margin (deg) 

44.1922 

34.4364 

Maximal tolerable delay (s) 

0.352 

0.237 


Fig. 2.6 shows the closed-loop step responses of these two systems with 0, 150 
and 300 ms added delay. With zero added delay, system I (upper figure) has slower 
responsiveness than system II because its bandwidth is lower. However, with a 150 ms 
delay, system I responds faster than system II, showing that a dynamic system with 
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higher bandwidth tends to be impacted more by the same amount of delay. With a delay 
of 300 ms, system I becomes more oscillatory but is still stable, whereas system II is no 
longer stable. Further analysis shows that although the two systems have the same phase 
margin, system I can tolerate much longer maximal time delay than system II (344 ms vs. 
238 ms). This can be interpreted using Eq. (2.12). Usually, the system with higher 
bandwidth tends to have higher crossover frequency {a > c ). It follows from Eq. (2.12) that 
the system with the higher crossover frequency suffers a larger phase lag <p d with the 
same time delay. And because^ = <p PM jo) , given the same phase margin ( <p PM ), the system 
with lower crossover frequency tolerates longer delay. 


Step Responses of System I (BW=4.3238 rad/s) with Delay of 0, 200 & 300 ms 



Step Responses of System II (BW=4.4864 rad/s) with Delay of 0, 200 & 300 ms 



Fig. 2.6. Step responses of the two dynamic systems to different delay 
2.2. Sources of Transport Delay 

This section is a summary of the sources of transport delay in a flight simulator, 
and some proposed methods to reduce the delay from each source. 
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Fig. 2.7 is a view of part of a flight simulator cockpit. The out-of-the-window 
display, i.e., the sea, the aircraft carrier and the lead helicopter, provides the visual cue to 
the pilot; the gauges, dials and photodiodes on the panel in front of the pilot provide the 
instrument cues, such as the airspeed, altitude and attitude; and the dynamic seat and the 
motion of the cockpit (if any) provide the motion cues. These are the three primary 
cueing channels, which serve as feedback to the operator in response to his input. There 
may be other cues available, e.g., the sound system may indicate the aural environment 
surrounding the pilot. 



Fig. 2.7. Cockpit view of a flight simulator 


When the pilot issues a command, the basic cues provide feedbacks to the 
operator, who continues to make control inputs based on these cues. The operator loop 
closes the man-machine system — the flight simulation. Because the simulator computers 
need time to generate the cues, they do not respond to the operator’s command 
instantaneously. The time it takes for the simulator to generate the basic cues is called 
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transport delay. Although the processing speed of computers has been improving, the 
increasing complexity of the simulator display negates the potiential to lessen the 
transport delay. As a result, transport delay still exists on state-of-the-art-flight simulators. 
This transport delay has a number of sources: 

1) Sampling delay — the time the simulation computer takes to sense a control input; 

2) Plant delay — the time the simulation computer uses to calculate the aircraft states; 

3) Data transfer delay — the time lapse from when the aircraft states are available to 
the point when the motion system, the visual system and/or the instrument system 
receives the signal; 

4) Cueing delay — the time it takes each individual system to respond to aircraft 
states. 

All these systems are likely to contribute transport delays, though these delays are 
not always present in concert, and more than likely, they will all be different from each 
other. 

2.2.1. Sampling Delay 

In a digital simulation, sampling delay arises from the fact that the simulator 
regularly senses the control input at the beginning of each frame, whereas the input 
arrives stochastically. For the zero-order hold (ZOFI*) system, the sampled value is held 
till the next sampling. Therefore, an input change in the middle of a frame is either lost or 


Though the accuracy of the first-order hold is higher than that of the ZOH, ZOH is used in the flight 
simulator because the first-order hold introduces an extra frame of delay. 
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delayed until the beginning of the next frame. The amount of delay varies in each frame, 
resulting in a sampling uncertainty. Fig. 2.8 illustrates the cause of sampling uncertainty. 


Input 



Fig. 2.8. Illustration of sampling delay 

The input change at point A coincides with the first sampling by chance, and thus 
the sampling delay is zero, representing the best-case scenario; the input at point B is not 
sensed by the computer until the second sampling, and the sampling delay is half a frame; 
point C is right after the 3 ld sampling, and it will be sensed in the 4 th sampling with a 
delay of almost a full frame, representing the worst-case scenario; for input at point D, 
the sampling delay is just a small part of a frame. With the exception of the best-case 
scenario, the sampled value represents the past input. Because the input may arrive at any 
time point during a frame, some sampling uncertainty results. As the number of frames 
becomes sufficiently large, the statistical average of the sampling delay is a half frame. 

Although the sampling delay can be reduced by using shorter frames, the frame 
length is driven primarily by the aerodynamics computation, rather than the sampling. 
Another way to reduce the average sampling delay is the multi-rate sampling proposed by 
R. M. Howe , but this results in higher hardware and software costs. 
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2.2.2. Plant Delay 


The plant in a flight simulator refers to the process of sampling the pilot control 
inputs, and then calculating the updated aircraft state using a specific aerodynamic model. 
The model first calculates the aircraft accelerations in all six degrees of freedom (DOF) 
in the body frame, then integrates the accelerations to get the velocities and 
displacements, and finally transforms the updated aircraft states (displacements, 
velocities and accelerations) to other frames, such as the world frame and the topodetic 
frame, for use by the individual cueing systems. Usually the process of calculating the 
accelerations, using the non-linear acceleration functions, is computationally intensive. 



Fig. 2.9. Bode diagrams of several numerical integration algorithms 

Although the two numerical integration processes themselves take little time to 
perform, they introduce phase lag and latency. Different integration methods have 
different characteristics in introducing phase lag, which is a function of the working 
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frequencies. Fig. 2.9 presents a comparison of some numerical integration methods in 
terms of both phase lag and gain distortion 4 . 

When two numerical integration methods are carried out one after another, they 
may cause latency due to the dependency of the latter on the outputs of the former for 
several iterations. For example, two 2 nd -order Adams-Bashforth integrations cause an 
extra frame of latency. If the second integration (from velocity to displacement) is 
replaced with the trapezoidal integration, this extra frame of latency could be removed. 
Therefore, the aerodynamics processing delay can be reduced by judiciously choosing the 
integration combination, and by rearranging the order of computation within each 
integration frame. 


2.2.3. Data Transfer Delay 

Data transfer time is the time it takes for the visual system (or the motion or 
instrument systems) to get the output signals computed by the simulation computer 5 . If 
the update rate of the simulation computer and the visual system are not the same and the 
latter is not an integer multiple of the former, communication asynchrony occurs which 
results in delay. Asynchronous delay occurs when the vehicle states are available from 
the host computer but the image processor is not ready to receive them until the 
beginning of the next frame. According to Richard McFarland 6 , the asynchronous delay 
between an image processor and a simulation computer appears to be random, but it is 
actually periodic. The periodicity can be demonstrated with the time lines shown in Fig. 
2.10. The periodic component of the asynchronous delay may be given as 


t d {m) = T 




(2.13) 
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where m is the iteration index of the visual computer, T is the update period of the host 

computer, and ^ is the minimal fraction equal to the ratio of (T is the update period 

of the visual computer), or M and N must be relatively prime. For the above example 
t d (1) = 30, t d (2) = 20, t d (3) = 10, t d (4) = 0, t d (5) = 30, • • • 

And the average asynchronous delay is 

( 2 - 14 ) 

iv m=] 

For the above example, 

T d = -^-(30 + 20 + 10 + 0) = ! 7. 5(/ns) 


Mainframe time line 

1 2 3 4 5 

Q Q O Q O 


Visual system time line 
E F 
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t d (3) t d (4) 
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Fig. 2.10. Illustration of periodical asynchronous time delay 

If the transfers are asynchronous, the data transfer delay affects the sampling 
delay. As long as the transfer time is less than the sampling interval (i. e., the frame 
length), transfer time may be considered the same as processing time 7 . To reduce or even 
eliminate the asynchronous time delay, the update rates of the subsystems must be equal 
to each other, or at least the update rate of the down stream subsystems should be an 
integer multiple of that of the upstream one. 
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2.2.4. Cueing Delay 

When the aircraft states are available from the simulation computer, the three 
basic visual, motion and instrument systems use them to generate the individual cues in 
the hardware systems, as depicted in Fig. 2.1 1. Each of the three cueing systems requires 
a different amount of time to generate and execute the hardware commands, which leads 
to cueing mismatches. 

A. Instrument Delay 


The quantities to be displayed on the cockpit instruments, such as the position and 


orientation of the aircraft, are available from the simulation computer, and hence there 


are very few time-consuming computations in the instrument system. The calculations 


necessary to transform the variables into the positions of dials of the indicators, and the 


digital-to-analog conversion contribute no more than one frame of transport delay. 



B. Motion Delay 


The purpose of the motion system of a flight simulator is to impart motion cues to 
the pilot, which are very similar to those he would experience in an actual aircraft. 
Because of the physical limits of any ground based simulator motion system, (i.e., it can 
not reproduce the large motions of a real aircraft), the simulator motion commands 
(attitude and rate) must be calculated from the aircraft states. This process is the core of 
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the motion cueing algorithm, which makes use of an operator motion sensation model 
and control theory. Various combinations of operator motion models and control theories 
lead to different types of motion cueing algorithms. The more advanced motion cueing 
algorithms may need to solve a Riccati equation that is computationally intensive. Telban 
and Cardullo have developed a neuro-computing approach to solve the Riccati equation 
in real time. After the motion base trajectory is computed, the simulator motion 
commands are then transformed from the degree-of-freedom space to actuator space 


(kinematic transformation), and the desired commands to the six actuators are generated. 

o 

A block diagram illustrating the function of the motion system is given in Fig. 2.12 . 



Simulator 

Motion 


Fig. 2.12. Block diagram of a motion cueing system (courtesy of Telban) 

The motion cueing algorithm and the kinematic transformation may take as long 
as one frame to complete, and the D/A conversion takes half a frame on average. 
Therefore, the motion system can cause as much as one and half frames of transport delay, 
although, usually the motion cueing algorithm is executed in the same frame as the flight 
dynamics. 
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C. Visual Delay 


The visual system consists of four parts in series: the front-end computer (FEC), 
the geometry processor (GP), the display processor (DP) and the display. The FEC takes 
care of the scene management and retrieving the image from the database; the geometry 
processor calculates the geometric transformations and the perspective transformation; 
and the display processor pixelates the data. Each of these three may consume one full 
frame. The fourth part (display) consumes a half frame on average, with the worst case 
being a full frame. Fig. 2.13 illustrates two scenarios of how these four steps in the 
pipeline may be implemented. In the worst case they work serially without overlap, 
causing 4 frames of transport delay; in the best case, parallel computation is employed 
between the FEC, GP and the DP, reducing the visual delay by one frame. 
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Fig. 2.13. Two scenarios for delay in the visual system 
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D. Cueing Mismatches 


The previous analysis shows that the instrument system causes the shortest 
transport delay, the motion system the medium, and visual system the longest. The 
discrepancy in the delay between the cueing systems is referred to as cueing mismatch. 
Fig. 2.14 illustrates the cueing mismatch between the motion and the visual systems. 
Frame 1 is the initial state; in frame 2, the pilot initiates a control input, but neither the 
motion nor the display changes; in frame 3, the response of the dynamic seat begins, but 
the display still remains unchanged; and finally, in the 4 th frame, the response of visual 
display begins. Though the research is not clear on the maximum mismatch a human can 
tolerate in a flight simulator, the cueing mismatch is believed to be the main cause of 
simulator sickness. Therefore, steps to sufficiently reduce the mismatch are necessary. 
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Fig. 2.14. Cueing mismatch between the motion and the visual systems 
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2.2.5. Summary 

Therefore, a typical flight simulator may incur, on average, 1.5 frames of 
transport delay in the simulation computer, and 3.5 frames of delay in the visual system. 
Assuming the two subsystems have the same frame length, and assuming an average 0.5 
frames of data transfer delay, the average total transport delay is then 5.5 frames, and 
may vary from 4 to 7 frames. 

If the transfers are asynchronous, the data transfer delay affects the sampling 
delay. As long as the transfer time is less than the sampling interval (i. e., the frame 
length), transfer time may be considered to be the same as processing time. Although the 
simulator time delay consists of several components from different subsystems, they 
make no difference to the operator, who only feels the total effect. 

2.3. Effects of Transport Delay 

It follows from the theoretical analyses on time delay in section 2.1 that time 
delay introduces oscillation to the closed-loop control system. In a man-machine system 
such as a flight simulator, this causes pilot induced oscillation (PIO). Increased 
oscillation means more and larger swings in response, which also mean larger variance of 
system output and degraded simulation performance. Conversely, the system output is the 
feedback on which the operator action is based. Increased output variance makes the 
operator compensate for the error, and this compensation causes the operator workload to 
go up. And for this reason, the operator’s perception of handling quality is degraded. 
Therefore, three effects of time delay on the flight simulation are expected: 

1) The system performance is degraded; 

2) The operator’s workload is increased; 
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3) The operator’s perception of handling quality is decreased. 

The three effects of transport delay are usually evaluated with three different 
metrics: namely, the Root Mean Square Error (RMSE) for the system perfonnance, the 
Power Spectral Density (PSD) for the pilot workload, and the Cooper-Harper Rating 
(CHR) for the handling quality. 

The RMSE is a measure of variance of the simulator output, and is used as an 
objective metric of the system performance. Larger RMSE means poorer perfonnance. 
The RSME is usually applied to simulation tracking tasks, such as compensatory tracking, 
pursuit tracking, or a gliding task. 

The power spectrum or power spectral density (PSD) of a digital signal is the 
discrete Fourier transfonn (DFT) of the autocorrelation sequence of a signal. In 
mathematical terms, the PSD is proportional to the square of the magnitude of the process. 
Hence, it is closely related to the energy of a signal as a function of the frequency, or it 
represents the frequency distribution of energy of a signal. The PSD of the pilot control 
input deflection (the pitch stick, roll stick, rudder pedal or throttle) is related to the energy 
the pilot put into the simulator, often referred to as the workload. It is also an objective 
metric. 

The Cooper-Harper rating scale 9 includes ten values, listed in Table 2.2, in 
ascending order corresponding to decreasing handling qualities, i.e., a scale value of one 
represents the best handling quality, and a value of 10 represents the worst deficiency. 
Ratings of 1, 2 and 3 fall in level I, ratings of 4, 5 and 6 are in level II, and ratings of 7, 8, 
9 and 10 are in level III. Jumping from one level to another is considered to be a 
significant change of handling quality. A CHR is a quasi-objective evaluation because the 
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operator is required to make a series of decisions concerning the difficulty of controlling 
the aircraft. The Cooper-Harper rating is an evaluation of the overall system, rather than a 
single factor such as responsiveness; hence two simulations with close CHRs may have 
quite different dynamic responses. 

In summary, the RMSE is an objective metric of the simulator output, the CHR is 
a quasi-objective metric from the operator’s cognitive judgment, and the PSD is an 
objective metric of the simulator input (or the operator output). 

Many researchers have extensively explored the effects of time delay on a flight 
simulation using these three metrics. A study of the available literature is summarized as 
follows. 

Table 2.2. The Cooper-Harper scale 



By analyzing the altitude errors in different delay conditions, Crane 10 reported 
that pilot dynamic response and system performance change as the pilot attempts to 
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compensate for the decrease in system stability caused by the transport delay. More 
importantly, the changes in pilot dynamic response and the system perfonnance can bias 
the results of a simulation by influencing the pilot’s rating of handling qualities of the 
simulated aircraft. 

Bailey, et al 1 1 , conducted some experiments on a simulator including four 
simplified aircraft models — fighter, small cargo, medium cargo and large cargo, with 
pilots flying some demanding tasks in both ground based and in-flight modes. 0-240 ms 
artificial delay was added over the 100 ms baseline delay. Bailey noted: 

1) Except for the fighter in the in-flight mode, the time delay significantly increases 
the RMSE, or degrades the performance. 

2) The time delay causes the regression line of the Cooper-Harper rating to cross the 
Level I and Level II boundary, indicating significant handling quality degradation. 

3) The difference in the effects of time delay among the four aircraft models is small. 

4) In the in-flight mode, the negative effects of time delay are not as serious as in the 
ground-based mode. Because the in-flight mode could be thought of including the 
motion cue, it suggests that the motion cue makes the operator less sensitive to the 
time delay. 

Similar simulation experiments were conducted by Riccio, et, al , in which 0-350 
ms time delay was added in addition to a baseline delay of 50 ms, with the pilot 
controlling a fighter or a large transport aircraft to maintain a constant heading and 100 ft 
altitude over flat terrain, and with turbulence. In this experiment only the RMSE 
evaluation was carried out. Both heading RMSE and altitude RMSE increased 
significantly with the time delay. There was no statistically significant difference between 
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the two aircraft models, however, the RMSE curve of the transport appeared to have a 
larger curvature. 

Middendorf, et al , on the other hand, used the power spectral analysis to 
investigate the effects of the simulator time delay on flight control activities. The 
simulation was conducted on a fix-base (no motion cue) simulator. Time delays of 90, 
200 and 300 ms were inserted in addition to the 90 ms baseline delay. The lateral stick 
movement was recorded for a spectral analysis while the pilot was completing an offset 
approach. For this task, the aircraft was initially lined up with one runway, then 
transitions to an adjacent parallel runway within a certain distance. They reported the 
following findings: 

1) There was a peak in the power spectrum at approximately 0.08 Hz, which is, as 
they stated, related to the inverse of the period of the sidestep maneuver. 

2) There was another peak at approximately 0.25 Hz, appearing to be a direct result 
of the maneuver itself. 

3) The reduced damping ratio resulting from the time delay made the closed-loop 
system less stable at high delay conditions. 

Cardullo conducted two phases of simulator experiments with his graduate 
students, Telban and Guo 14 , and some colleagues at the NASA Langley Research Center. 
The experiments were conducted on the Visual Motion Simulator, and they were 
designed to test novel motion cueing algorithms as well as the effects of time delay and 
the McFarland compensation algorithm. Artificial delays of 0, 50, 100, and 200 ms were 
added to the simulation. These delays, combined with the three motion cueing algorithms 
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— the adaptive, optimal and nonlinear algorithms, generated 16 different conditions. Each 
pilot was assigned three tasks — a straight-in approach, an offset approach and a takeoff. 
Half of the takeoff test runs included an engine failure. Three pilots took part in the first 
phase of the experiment (also referred to as the preliminary experiment), and eleven took 
part in the second phase. Cooper-Harper ratings were logged only in the preliminary 
experiments. In both phases of the experiment, the time domain data of the pilot’s control 
deflections and some other simulator state variables were recorded for power spectrum 
analysis. Although most of the results are outside the scope of this section, the results on 
the effects of time delay and compensation are summarized here. 

1) The time delay slightly increased the overall integrated PSD of the roll control, 
pitch control and rudder pedal control, and significantly increased the integrated 
PSD of the roll and pitch controls in the frequency range of [0.17 0.4] Hz. 

2) The time delay moved the highest peak of the PSD plot to higher frequency area. 
The effects of the engine failure in the takeoff maneuver dominated the total PSD 
of the pilot control inputs, obscuring the effects of the time delay on PSD. 

3) The Cooper-Harper ratings did not increase with time delays of up to 100 ms, in 
the two approach maneuvers. 

Transport delay may also cause simulation sickness. The tenn “Simulation 
Sickness” is usually reserved for situations that are nauseogenic in the simulator but not 
in the corresponding aircraft. Simulation Sickness shares many symptoms, such as 
drowsiness, dizziness nausea, vomiting, etc, with the motion sickness that occurrs in the 
real systems. It is difficult to conduct simulator sickness research because it is difficult to 
objectively evaluate it. Literature research suggests that cue mismatch is the primary 
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cause of simulator sickness. Duh, et al, recently proposed a visual-vestibular crossover 
frequency concept, and hypothesized that conflicting visual and inertial self-motion cues 
at the frequency of maximum crossover would be more likely to evoke simulator sic kn ess 
than conflicting cues at a high frequency. This hypothesis is supported by experiments. 
Sobiski states, “Temporal cue mismatch can contribute to a malaise known as ‘simulator 
sickness’ and it may be due to the frantic efforts of the brain to resolve two conflicting 
sources of motion related information.” Because the cue mismatch is caused by the 
difference in time delay, and simulator sickness is induced by the cue mismatch, it may 
be inferred that simulator sickness is indirectly caused by the difference in time delay. It 
seems logical that time delay may cause simulator sickness. Long time exposure to PIO 
caused by time delay is expected to make the operator feel tired and to induce some 
instantaneous or simulator sickness. However, to date, the direct relationship between 
time delay and simulator sickness is not completely understood, probably because there 
has been insufficient research on this issue. Zaychik 15 reported that Nelson et al, 
investigated the influence of time delay, time on task and task complexity on subjective 
ratings of simulator sickness. It was shown that the subjective ratings varied directly with 
duration of task. The delay variance affected the performance and workload of a subject 
but had no effect on the simulator sic kn ess questionnaire (SSQ) ratings. Zaychik suggests 
that this interesting observation can be attributed to the subjective nature of the simulator 
sickness assessing technique used by the researchers. 

Likewise, Uliano, et al 16 , reported that visual lag had no effect on illness in their 
experiments. Although longer lags were somewhat disruptive of perfonnance, there was 
no evidence that they contributed to illness. These results are surprising because the range 
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of asynchrony studied was fairly large (126 ± 17 ms to 215 ± 70 ms). Based on these 
results, Uliano et al suggested that within the range of operationally useful simulators, 
visual asynchrony does not appreciably contribute to simulator sickness. Lags 
approaching 300 ms in flight simulators become too unrealistic and/or too difficult to fly. 
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3. Current Techniques of Compensating Transport Delay 

As stated in Chapter 2 (section 2.3), transport delay in a flight simulator reduces 
the system phase margin, and can cause pilot-induced oscillation (PIO), which in turn 
leads to poorer system performance and handling quality, higher control workload, and 
even simulator sickness. Transport delay may be reduced, by rearranging the order of the 
calculations, choosing more efficient algorithms, increasing the sampling rate, or 
synchronizing the communications. If the delay is still exceeds the tolerable level, for 
maintaining desirable simulator perfonnance, algorithms to compensate for the delay 
may be employed. The purpose of compensation of time delay, from a time-domain 
perspective, is to provide prediction to counteract the system latency; and in terms of the 
frequency-domain, is to create phase lead to offset the reduction in system phase margin. 

When transport delay exists in a visual system, the image that is displayed is 
actually the delayed image representing a past aircraft state. Since one cannot generate an 
undelayed image from a delayed one, compensation must be applied to the aircraft state. 
Therefore, the idea of compensation is to predict the future aircraft state using the 
currently available state information. Therefore, an image, based on the predicted aircraft 
states, can be used to offset the transport delay in the visual system. This idea is 
illustrated in Fig. 1.4, where the dashed line in the small plot to the right of the predictor 
block, illustrates the prediction. 

The dashed predicted curve in Fig. 1.4 is an ideal prediction, i.e., it has exactly the 
same shape as the original solid one, with only a pure phase shift. Unfortunately, an ideal 
predictor does not exist in the real world. A practical predictor generally introduces error, 
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and the longer the time delay, the greater the error. Mathematical analysis illustrates that 
pure time delay brings phase lag to the system without changing the magnitude (Section 
2.1). Therefore, a good predictor must satisfy two basic criteria: 

1) It must be able to provide sufficient phase lead to offset the phase lag caused by 
the time delay. 

2) It must introduce minimum gain distortion. 

In addition, the compensation must be simple enough that it does not introduce 
extra time delay due to computation. Therefore, a third criterion is: 

3) The computation workload of the predictor must be minimal. 

To date, many compensation techniques have been developed to mitigate the 
transport delay in the flight simulator. The lead/lag filter, the McFarland compensator and 
the Sobiski/Cardullo state space predictor are the three most prominent current 
techniques, and will be considered here. 


3.1. Lead/lag Filter 


The transfer function of the lead/lag filter is given by 
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(3.2) 


where Y p (s) and Y(s) are the Laplace transforms of the predicted aircraft state and the 
undelayed aircraft state, respectively; and (O n and co d are the two corner frequencies of 
the filter. If co n < co d , the Bode asymptotes of both the magnitude and the phase of the 


filter given by Eq. (3.1) are illustrated in Fig.3.1. 
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Fig. 3.1. Bode asymptotes of the lead/lag filter 

The hump of phase asymptote in \a) n co d \ provides the phase lead to the system, 

and the maximal phase lead occurs at a frequency (O m called the medium frequency, 
which is the geometric mean of the two corner frequencies as given by 

a> m = (3-3) 

and the maximal phase lead is 

• , 1 -cc .. 

sin^ = (3.4) 

1 + a 

This phase lead is obtained at the expense of gain distortion because the 
magnitude \H {joA is not unity. And, since the high-frequency gains are increased, any 

system using phase lead compensation may be subjected to high-frequency noise 
problems. Nevertheless, as Crane 10 states, the resulting increase in system gain at 
frequencies greater than the crossover frequency ( (O c ) is not normally a problem, because 
the system amplitude ratio and the power of the input and disturbance signal usually 
decrease rapidly at frequencies greater than a> c . 

Designing a lead/lag filter involves choosing the gain k and the two comer 
frequencies co n and co d (or the two time constants r„ and z d ). The lead/lag filter had long 
been used in other control systems to provide phase modification before it was first 
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applied in a flight simulator. Using the classical frequency-domain design method can 
meet the design specifications accurately. However, because the classical method 
assumes linearity, and the aerodynamic models of a flight simulator are usually non- 
linear, the classical approach does not lend itself to the design of a lead/lag filter for 
compensating the transport delay in the flight simulator. 

Because the pilot crossover frequency region has been shown to be the most 
critical for pilot control, and for pilot ratings (Fig. 3.2 17 ) of the fidelity of a dynamic 
simulation, an ideal design would be to make the medium frequency co m , the frequency at 

which the maximal phase lead occurs, at the pilot crossover frequency co c . Unfortunately, 
co c is usually unknown. One approach would be to assume an estimated OJ , and let 
co m = &> c . Then calculate the maximal phase lead necessary to counteract the phase lag 
caused by the delay t d 

</> = t d G) c (3.5) 

Then calculate the ratio a according to Eq. (3.4), and finally determine the two time 
constants. Though simple, this design approach has two problems. First, the estimated 6t) c 
is not the real crossover frequency (O c , which varies in a simulation due to many factors. 
The further a> c (i.e., co m ) departs from oj , the less phase lead at oj compared to the 
maximal lead at co m . In other words, the phase compensation at the crossover frequency 
is not sufficient. Second, since 0 < or < 1 , it follows from Eq. (3.4), that (f) must be less 
than tt/ 2 . Then, according to Eq. (3.5), for a long delay t d , the estimated OJ must be 
very small, and hence it also diverges from co c , resulting in insufficient phase 
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compensation. This explains why the lead/lag filter cannot be used to compensate long 
time delay. 



Fig. 3.2. Pilot sensitivity envelops in the frequency domain 

Therefore, the design of a lead/lag filter to compensate the time delay in a flight 
simulator with unknown or nonlinear aerodynamic model cannot be exact, and therefore 
approximate methods are usually used. Two such approximate methods were introduced: 
one by Crane, and the other by Ricard and Harris . In both methods, the numerator time 
constant T n is set to be equal to the time delay t d , whereas the denominator time constant 

T d is chosen in a different manner. In Crane’s approach, z d is the solution of the 
following equation 


°/L=«. = tan ~ X ( a c T d ) = ( 3 - 6 ) 

The gain k can be chosen so that the filter gain is unity at the crossover frequency. Eq. 
(3.6) then becomes problematic because, for the widely-accepted value of co c of 2-3 Hz 
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and t d up to 800 ms, the term tarf 1 [co c t d ) is always less than 0 ) c t d . Hence the term 

tan ' (o\t a ) is negative. Any T d that makes tan ' (co c T d ) negative is either negative 

(meaningless) or greater than T n (in this case, the filter becomes a phase lag filter). As 

stated previously, the actual phase lead of the filter is usually less than the maximum, and 
one way to resolve this difficulty is to multiply the right side of Eq. (3.6) by a coefficient, 
thereby changing it to 

= tan ' (MJd)-*™- 1 ( 0) c T d ) = T]co c t d , (0 < 77 < 1) (3.7) 

where 77 is dependent upon the crossover frequency and the time delay. For example, for 
a crossover frequency of 2 rad/s and time delay up to 800 ms, 77=0.5 produces very good 
results. Table 3.1 gives the design results obtained using the modified Crane method for 
time delays of 200, 400 and 800 ms and an estimated crossover frequency of 2 rad/s. 

Table 3.1. The lead/lag filter coefficients designed with Crane’s method 


t(l (ms) 

C, 


a>m (rad/s) 

(deg) ; 

200 

0.2 

0.0912 

7.4044 

21.94 

400 

0.4 

0.1410 

4.2108 

28.60 

800 

0.8 

0.1080 

3.4021 

49.65 


In Ricard’s approach, T d is chosen such that the best pilot flight control 

performance is achieved. Using Ricard’s method, the two time constants for lead/lag 
filters used to compensate for delays of 200, 400 and 800 ms are listed in Table 3.2 
(obtained from Sobiski 19 ). Column 4 gives the frequencies at which these filters provide 
the maximal phase lead. The maximal phase lead is listed in column 5. Assuming the 
crossover frequency of the flight simulator system is 2 rad/s, the decrease in phase 
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margin caused by time delays of 200, 400 and 800 ms are 22.9°, 45.8° and 91.7°, 
respectively. These values are greater than the corresponding maximal phase lead of the 
lead/lag filter (column 5 in Table 3.1). This shows that none of these filters provide 
sufficient compensation. Although other design methods may yield better compensation, 
the improvement is limited. The lead/lag filter simply lacks the ability to provide large 
phase lead. 

Table 3.2. The lead/lag filter coefficients designed with Ricard’s method 


t d (ms) 

T n (S) 

T d (s) 

CO m (rad/s) 

<t> m ( de g) 

200 

0.2 

0.1859 

5.2 

2 

400 

0.4 

0.2105 

3.5 

18 

800 

0.8 

0.1695 

2.5 

40 


Reguarding the effects of the lead/lag compensation, Crane reported, “The 
compensation is effective; improvements in pilot performance and workload or HQR 
were observed. The delay compensation approach attempts to minimize changes in pilot- 
aircraft dynamics in the region of the crossover frequency.” 

Ricard and Harris reported that there were effects due to the presence of the 
lead/lag filter, and the error score indicated that varying z d was significant, but not 

significant for the crossover power. However, Ricard and Harris did not reveal whether 
the pilot performance, similar to that in the undelayed system, was achieved when using 
lead/lag compensation. 

To test the effectiveness of the lead/lag filter in compensating the transport delay, 
a lead/lag compensator was added to the control system shown in Fig. 2.3 in front of the 
added time delay. The new system block diagram is given as Fig. 3.3, where y is the 
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undelayed aircraft state, y p is the predicted aircraft state and y c is the delayed prediction, 

or the compensated state. The added time delay is 200, 400 or 800 ms, and for each delay 
case, both the lead/lag filters designed with the modified Crane method and the Ricard 
method (see Table 3.1 and Table 3.2) were applied. The frequency responses and step 
responses are given in Fig. 3.4 and Fig. 3.5. 



Fig. 3.3. Block diagram of a delayed control system with a compensator 


CO 

■o 

CD 

10 

TJ 

D 

’E 

0 

CD 

CD 

-10 


-20 


10 


Bode Diagrams of a Control System with Delayand Compensation 
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Fig. 3.4. Bode diagrams with delay (400 ms) and with/without lead/lag 
compensation of different gains 

Fig 3.4 shows the 400 ms delay compensation using the modified Crane design. 
When the filter gain is unity (dashed curve), it provides 20 “phase lead at 2 rad/s. But 
because of the gain distortion, the crossover frequency is shifted to 3 rad/s, corresponding 
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to a phase margin of -38.44°, which is even lower than the delayed system without 
compensation. To move the crossover frequency back to 2 rad/s, the pilot can decrease 
his gain, or the filter designer can decrease the lead/lag filter gain from unity to 0.833 
(dotted curve). The resulting lead/lag filter provides a phase lead of 12.43°, which is 
equivalent to a 108 ms delay, and much less than desired 400 ms. 


Step Responses: Undelayed, Delayed and Compensated 
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Fig. 3.5. Step responses of the undelayed, delayed and compensated systems 

From the step responses, a lead/lag filter designed using Crane’s method achieves 
better results than one designed with Ricard’s method in all delay cases. The results are 
greatly improved for a delay of 200 ms, in which Ricard’s filter provides little phase lead. 
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In all cases, a lead/lag filter cannot compensate the entire delay: no matter how long the 
delay is. The lead/lag filter is only beneficial when used to compensate for delays of less 
than approximately 100 ms. This should be considered the upper limit of the lead/lag 
filter. 

In short, although the lead/lag compensator does meet the third criterion of a good 
compensator, it does not meet the first two criteria.. 

The pulse transfer function and the difference equation of the lead/lag filter are 
given by Eq. (3.8) and (3.9), respectively 


Y p{z) _Z + C n 
Y(z) z + c d 

(3.8) 

c d y P {k) + y{k + l) + c n y{k) 

(3.9) 


As an Infinite Impulse Response (HR) filter (since c d ^ 0 ), the lead/lag 
compensator makes use of the previous prediction to calculate the current prediction, thus 
the error from one iteration is passed on to the next iteration, resulting in error 
accumulation. This error accumulation is the primary cause of the gain distortion. 

3.2. McFarland Filter 

20 

Richard McFarland developed a Finite Impulse Response (FIR) filter to avoid 
the weaknesses of the lead/lag filter that cause the prediction error to accumulate. Its 
pulse transfer function and difference equations are given by 

Y p (z) = Y(z)+(b a +b l z- l +b 1 z- 1 )v(z) (3.10) 

y, (k) = y(k)+b ( ,v(k) + b i v(k-i) + b,v(k-2) (3.11) 
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where y is the aircraft state to be predicted, v is the corresponding velocity, y p is the 
predicted aircraft state, k is the iteration index, and 7 (z), Y (z) and V (z) are the z- 
trans forms of y p , y and v , respectively. Clearly, the McFarland compensator is a 
special integrator making use of three consecutive steps of velocity. The three 
coefficients b 0 , b x and b 2 determine its ability to compensate the time delay. McFarland 
uses a method known as sinusoidal tuning to detennine these three coefficients, which are 
the solutions to three equations derived from the boundary conditions of the pass band 
[0 co 0 ] by assuming a sinusoidal input to the fdter. The pass band comes from the 
assumption that the pilot operates primarily within this frequency interval and the 
operation beyond co 0 , about 6-20 rad/s, is insignificant. At the zero frequency condition, 

the velocity is constant, and Eq. (3.11) becomes y (k) = y(k) + (b 0 +b x +b 2 )v(k ) . For 
an ideal prediction of y (k) that is t d ahead of y(k) , the relation 
y (k) = y(k) + v(k)t d holds. Comparing this with the previous equation, and the first 
equation from the boundary conditions is obtained 


b 0 +b x +b 2 =t d 

*d_ 

In addition, by substituting the relations Y p (z) = z T Y(z) and Y (z) 


(3.12) 


T\ 

n-z- 1 ^ 

2 

li +z_1 J 


y(z) 


(trapezoidal integration) into Eq. (3.10), the relationship becomes 

H*) 


r(z) 


: Z < — 


1 + Z 

f 


-1 


+ (6 0 +b x z 1 +b 2 z 2 ) 


(3.13) 
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Let Y(z)/V(z) be equal to the ideal integration 1 / (jco) at the other boundary frequency 


co 0 , and the magnitude and phase angle equalities give the remaining two equations. 
Combining these with Eq. (3.13), and the three coefficients are given as follows, where 
0 O = co 0 T and if/ 0 = co 0 t d . 


b n =- 


\j// Q + siny/ 0 (1-2 cos 9 0 ^sin 0 O + k, G 0 sin6 0 -cosi/Zq (1 — cos 6 {) ) (l + 2cos 6 0 ) 


2 co 0 sin 6 a (\ — cos 0 Q ) 


(3.14) 


2sin(6 , 0 + ^ 0 )-2^ 0 cos# 0 -0 O (l + cos^ 0 ) 
2<y 0 (1 -cos^o) 


(3.15) 





sin O 0 - cos i/f 0 ( l- cos 0 O ) 


2<y 0 sin 0 O (1 - cos G 0 ) 


(3.16) 


McFarland states that “the algorithm delivers high-fidelity, compensated CGI 
drive signals over the human-factors bandwidth, and can dramatically improve the pilot 
control for high gain tasks such as precision hovering and station keeping.” However, he 
does not show experimental results to substantiate it. A system similar to the one depicted 
in Fig. 3.3 was used to test the McFarland compensator. In this test, the lead/lag filter was 
replaced with the McFarland filter. The system was then transfonned into a discrete form 
before the compensator was applied because the McFarland filter was in discrete fonnat. 
The frequency responses and step responses are given in Fig. 3.6 and Fig. 3.7, 
respectively. 
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Bode Diagrams of Systems with No Delay, 0.2s Delay, McFarland Compensation 



Frequency, rad/s 


0 

-200 

O) 

CD 

■o 

a) -400 

CO 

03 

_C 

-600 

-800 


i 1 1 i ill 

1 1 1 i i i i i 1 i 

• 1 ' 



i i • i i i i i i i 

! !!!!!!! 

i i i i f 4 sYk i 


With No Delay 

i 

i 

1 

i 

- - 4 

•-* 

■ - • 4 

• 

• 

• 

'± 1 

- 

• With 0 . 2 s Delay 

1 fTTT m iT““ V 


McFarland Compensation 

1 1 1 1 1 1 1 1 1 

1 1 1 1 1 1 1 1 1 

1 1 1 1 l 1 1 1 1 1 


10 


10 


10 


Frequency, rad/s 


Fig. 3.6. Bode diagram of McFarland compensation for delay of 200 ms 

The Bode diagram shows that for time delay of 200 ms, the McFarland filter 
provides satisfactory phase compensation when the frequency is below 5 rad/s, but the 
phase lead is not sufficient in higher frequencies. The gain distortion is small when the 
frequency is within the pass band, but the gain distortion escalates when the frequency is 
beyond the pass band. For comparison, the compensation by the lead/lag filter designed 
with the modified Crane’s method is also plotted in Fig. 3.7. This plot shows that the 
McFarland compensator can provide more phase lead than the lead/lag filter, however, 
the McFarland filter still cannot provide 100 percent compensation. 
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Step Responses of the Systems of No Delay, Delayed by 0.2s and with Compensation 



Fig. 3.7. Undelayed, delayed by 0.2s and compensated responses 

The McFarland predictor can also be applied to simulation data. Fig. 3.8 shows a 
time history of roll angle and roll velocity recorded on the Visual Motion Simulator at the 
NASA Langley Research Center. By using the roll angle as y and the velocity as v , and 
with a prediction time of 192 ms, the McFarland compensator described by Eq. (3.11) 
produced the dashed curve shown in Fig. 3.9, this is the signal y in Fig. 3.3. The 

prediction is then delayed by the same amount for ease of comparision with the 
uncompensated roll angle. For ideal compensation, the delayed, compensated curve 
(dashed line) should be exactly the same as the undelayed uncompensated curve (solid 
line). However, this plot shows that the McFarland prediction has two problems: the 
actual prediction is only 176 ms, less than the design prediction by almost a frame (16.5 
ms), and more seriously, it causes very large spikes. 


50 


Recorded Roll Angle, Velocity from a Simulation in the VMS 
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Fig. 3.8. Roll angle and velocity of a real simulation 

The gain distortion and insufficient phase lead at high frequencies, demonstrated 
here, were verified by the perfonnance of the pilots in the simulations conducted in the 
VMS at the NASA Langley Research Center in 2002. Fig. 3.10 shows the power spectral 
densities of the roll control averaged across the pilots for both delayed and compensated 
(McFarland compensation) cases in an offset landing maneuver with 200 ms artificial 
delay added to the visual system. The McFarland compensator significantly reduces the 
PSD in the frequency range of [0.5 2.3] rad/s, whereas it increases the PSD in the range 
of [2.3 4.4] rad/s. In addition, the Cooper-Harper ratings indicate that for long added 
delay (200 ms), the pilot handling qualities tend to be degraded by the McFarland 
compensation. 
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McFarland Compensation to Simulation Roll Angle 



Fig. 3.9. Prediction by the McFarland filter of the simulation data 
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Fig. 3.10. PSD of the roll stick with and without McFarland compensation 

Apparently, the biggest problem with the McFarland filter is the annoying spikes 
it causes. The longer the delay is, the larger the magnitude of the spikes. The McFarland 
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algorithm is a special integrator or an extrapolator. The simplest extrapolator to provide a 
prediction of t d is 

y P ( k )=y( k ) +t d v ( k ) ( 3 - 17 ) 

If the velocity changes slowly (low frequency), this works well, but if the velocity 
changes quickly, the prediction introduces error because the velocity may be quite 
different t d later. The McFarland filter is superior to the extrapolator given by Eq. (3.17) 

because it uses three consecutive steps of velocities that can extrapolate the future 
velocity better than a single velocity. For moderate frequencies (around 1 Hz), let the 
average of these three velocities be v , then (3.11) reduces to y p (k) = y (k) + t d v (k ) , 

which is similar to Eq. (3.17), but the average velocity is used, resulting in a better 
prediction. Therefore, its working frequency range is wider than the pure extrapolator. 
However, if the velocity changes abruptly (the frequency is even higher), such as the 
velocity plotted in Fig. 3.8 near time = 60 seconds, spikes occur. No matter what the 
delay is, the three coefficients b 0 , b x and b 2 are always positive, negative and positive, 

respectively, and the absolute value of b x is always the largest, e.g., for t d =0.2s, 
b 0 =2.9979, b x =-5.5 197 and b 2 =2.72 19. The absolute values of the coefficients are at 
least 10 times greater than t d . If the velocity changes by more than 10 percent in several 

iterations, spikes are likely to occur. Table 3.3 gives an example, where the two spikes 
are highlighted. From this example, it can be shown that the spikes from the McFarland 
compensation are caused by: first, the three coefficients change sign alternatively, and 
second, they do not change value in relation to changing simulation conditions. A better 
choice for the three coefficients would be b 0 =b x =b 2 =t d /3, however, if the coefficients 


53 



could be made to change with the simulation conditions, the prediction could be even 
better. 


Table 3.3. Several iterations of McFarland prediction with spikes 


'( S ) 

v 0 (rad/s) 

Vj (rad/s) 

v 2 (rad/s) 

- Vp (rad) 

59.776 

0.1000 

0.1015 

0.1029 

0.0115 

59.792 

0.0877 

0.1000 

0.1015 

- 0.02 

59.808 

0.0754 

0.0877 

0.1000 

0.0085 

59.824 

0.0630 

0.0754 

0.0877 

0.0072 

59.840 

0.0507 

0.0630 

0.0754 

0.0056 

59.856 

0.0504 

0.0507 

0.0630 

0.0402 

59.872 

0.0502 

0.0504 

0.0507 

0.0080 


3.3. Sobiski/Cardullo Filter 

2 1 

In 1987, Sobiski and Cardullo proposed a state space predictor for compensating 
the transport delay. It is based on the equation 

x(t + t rf ) = e Md x(t) + 1 J e A ^ _r) Bw(t + r)dr (3.18) 

which is derived from the solution of the state space differential equation x = Ax + B u . 
This equation shows that the predicted state vector x(t + t d ) may be calculated from the 

current state vector x(t) provided that the future input u is known between t and t + t d . 

Unfortunately, this is an obviously impossible condition with stochastic operator’s 
control input u . Therefore, Sobiski made some assumptions about the form that the input 
might take, i.e., piece-wise constant, sinusoidal, exponentially decaying, etc, so that the 
future input may be approximated by the current input. Then the prediction is given by 
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x(t + t d ) = ^e x ‘ d ~\x(t)+ j ( ‘ e A<t '' t] dz B u{t) (3.19) 

By denoting 

<& = e A,d (3.20) 

and 

¥ = e A(,d ~ T) dr (3.21) 

Eq. (3.19) is simplified to 

x ;; =<Dx + ’FBw (3.22) 

3> and T are called the state transition matrix and the convolution integral matrix, 
respectively. Assuming the matrix D is zero (this is true for most control systems), then 
the predicted output is calculated by 

y p = (CO)x + (C*FB)m (3.23) 

Because C'FB 0 , the matrix D of the compensated system is usually not zero even 
though the matrix D of the undelayed system is zero. The problem caused by the non- 
zero D matrix of the compensated system will be addressed later in this section. Directly 
from Eq. (3.19), the structure of the Sobiski/Cardullo filter is illustrated in Fig. 3.11. The 
Sobiski/Cardullo state space predictive filter is an original approach for compensating the 
time delay. Theoretically it can compensate longer delay than the McFarland 
compensator because it uses more system information, i.e., the full state vector, though it 
requires complicated calculations. 

The state space filter may be expressed in discrete form. Let the discrete state 
equation be x(£ + l) = Gx(k) + Hu (k) , where G = e AT and H = f£ e AT dr\j$, where T 


55 



is the sampling period. The discrete state space filter that predicts / iterations in the 
future is then given by Eq. (3.24) 



Fig. 3.11. Sobiski/Cardullo compensator 


x(k + /) = [G / ]x(^)+ H u i k ) (3.24) 

L;=o 

In deriving Eq. (3.24), the same assumptions used by Sobiski reguarding 
continuous system input u are also held, so that the future inputs u(k + j) can be 

approximated by the current input u (k) . Define 

®,,(/) = G' (3.25) 

and 

'*’.(')= ( 3 - 26 ) 

j = o 

These are the discrete state transition matrix and convolution integral matrix, respectively. 
The drawback of the discrete state space filter (Eq. (3.24)) is that it can only predict delay 
in integer multiples of the frame time T . If the time delay t d cannot be divided exactly 

by T , say t d = IT + r, (0 < r < T) , the exact amount of prediction may be calculated with 
the interpolation method 

x((k + /)T + r) = [e Ar ]x((k + /)r)+ f'e^dA B u(k) (3.27) 
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Comparing Eq. (3.17) and (3.19), one can see that the Sobiski/Cardullo algorithm 
is a type of extrapolator in state space fonnat. The primary difference is that the 
Sobiski/Cardullo filter uses more system information, including the control input 
information. In the continuous state space predictor, more derivative or integration 
information is used. In the discrete state space predictor, more past information is used. 
Inclusion of the control input is a distinguishing characteristic of the Sobiski/Cardullo 
filter when compared to the lead/lag or the McFarland filter. 

To test the Sobiski/Cardullo compensator, the same system as depicted in Fig. 3.3 
is used, except the lead/lag filter is replaced with the Sobiski/Cardullo filter. Because it is 
in the state space format, the system is redrawn as in Fig. 3.12. The pilot model and the 
aircraft model are cascaded together to form the matrices A , B , C and D , and to 
calculate d> and T . As was done for the other two filters, 200, 400 or 800 ms delay and 
compensation are tested, and the frequency responses and step responses are shown in 
Fig. 3.13 and Fig. 3.14. 



Fig. 3.12. Block diagram of a control system with a state space compensator 


The Bode diagrams show that the Sobiski/Cardullo filter can compensate a 200 
ms delay with unnoticeable gain distortion and phase error up to 30 rad/s. In 
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compensating a 400 ms delay, the gain distortion is very slight, and the phase error is 
unnoticeable for frequencies lower than 10 rad/s, but the compensation displays slight 
phase lead that increases with frequency above 10 rad/s. In the 800 ms delay 
compensation case, gain distortion is significant in high frequency area, and the phase 
lead is insufficient for frequencies above 5 rad/s, causing the system to be unstable. For 
this reason, the step response of the 800 ms delay case is not presented. 


Bode Diagrams of the Undelayed System and the Compensated Systems (Sobiski/Cardullo) 
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Fig. 3.13. Bode Diagrams of the Compensated Systems (Sobiski/Cardullo) 

The step responses shown in Fig. 3.14 verify the frequency response analysis. The 
“compensated” response (thick dashed curve) refers to y c in Fig. 3.12, and the result 

obtained by delaying the prediction by t d is illustrated by the dashed dot curve. For good 

compensation, the compensated response must be as close to the response of the 
undelayed system (solid curve) as possible. The Sobiski/Cardullo filter achieves perfect 
compensation for 200 ms delay in this system, and for delay of 400 ms, the compensation 
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error is insignificant. Since these are all closed-loop step responses, they illustrate the 
change in responsiveness caused by the prediction and delay. Fig. 3.15 directly compares 
the lead/lag filter, the McFarland filter and the Sobiski/Cardullo filter to step responses. 
The superiority of the Sobiski/Cardullo filter is obvious in terms of both gain distortion 
and phase compensation. 


Step Responses of Undelayed System and Compensated System with S/C 




Fig. 3.14. Step responses of the Compensated Systems (Sobiski/Cardullo) 

Flowever, the Sobiski/Cardullo filter has several problems. First, since the state 
space equations are derived from linear differential equations, and they are constant in 
Sobiski’s implementation, the filter can only be applied to linear time-invariant (LTI) 
system. As previously stated, the aerodynamics of a flight simulator are usually nonlinear 
and time variant, and thus the matrices A , O and *P are not available. This is the 
primary reason the Sobiski/Cardullo filter has remained in laboratory use since its advent 
almost 20 years ago. Second, Sobiski’s implementation has some limitations. In his 
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implementation, the Pade approximation of the time delay is added before the 
compensator, as is shown in Fig. 3.11., (where Pade approximation of the time delay is 
necessary and the time delay model is cascaded with the operator model and the aircraft 
model so that the matrices O and ¥ carry the infonnation of the time delay) in contrast, 
the author’s approach is shown in Fig. 3.12. In Sobiski’s system, the delayed variable to 
be compensated must have the same dimension as the undelayed one. In a flight 
simulator, the transport delay appears to be the delayed image displayed on the screen. 
Because the image and the corresponding aircraft state do not have the same dimension, 
the real transport delay is not linear. In the author’s implementation (shown in Fig. 3.12) 
the delay occurs after the prediction which more closely represents the compensation in a 
real flight simulator. In a real system, the “Time Delay” may be a pure transport delay in 
the visual system, and may not necessarily match the Pade approximation of the time 
delay. Therefore, Sobiski’s implementation has more theoretical significance than 
practical usefulness. Third, the assumptions for approximating the future control input 
with the current one, such as piece-wise constant, sinusoidal, exponentially decaying, etc, 
do not apply to the real simulation conditions. As given by Eq. (3.23), the matrix D for 
the compensated system is not zero, and hence the high frequency components of the 
stochastic pilot control input u will be added to the compensation y . Finally, the state 

space filter involves a large computation burden because of the matrix operations, 
especially the power series operations in the calculations of the state transition matrix O 
and the convolution integral matrix ¥ . Because of these limitations it would be desirable 
to develop a more practical state space predictor that can compensate longer delay than 
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Roll Angle, deg 


currently available time delay compensators, and simplify the algorithm to provide 
minimal computation cost. 


Step Responses of the Systems of No Delay, Delayed by 0.2s and with Compensations 



Fig. 3.15. Comparison of the Three Prominent Compensators 
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4. Novel Approaches to Compensation of Time Delay 


As stated in Chapter 3 (section 3.2), there are numerous problems with the 
currently available compensators. The McFarland predictor offers insufficient phase lead 
and large gain distortion at high frequencies due to the constant coefficients that do not 
change with the simulation conditions. The Sobiski/Cardullo filter requires a linear time- 
invariant system model, whereas the aerodynamic models used in modern flight 
simulators are usually non-linear, time -variant and are frequently not readily available. 
Another major problem with the Sobiski/Cardullo filter is its computational burden. This 
chapter presents some new techniques to address these problems. First, a simple spike- 
reduction algorithm for the McFarland filter is introduced. Next, this chapter discusses 
using least squares methods in both frequency and time domains to design a three- 
velocity predictor. After that, the author will present two novel compensation algorithms 
that have been developed — an adaptive predictor which uses a Kalman filter and a state 
space predictor which uses a linear model of the aircraft dynamics to predict future states. 
The well-known Kalman filter technique is used in a unique manner so that the predictor 
can accurately provide the desired amount of prediction. From five different 
implementations of the Kalman estimator, the best option was selected based on the 
results of theoretical analyses. The state space filter with a linear reference model is the 
first practical model referenced state space predictor applied in a flight simulator to 
compensate the time delay. Several currently available linear aircraft dynamic models 
were tested, and the one that achieved the best prediction, based on the offline tests (see 
Chapter 5, pi 06), was chosen. The relationship between the reference model and the 
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quality of prediction was also investigated. By simplifying the state predictor to an 
ordinary predictor in a transfer function format, the computational workload is reduced 
significantly. The two new compensators are described below. 

4.1. Reduction of Spikes in the McFarland Compensator 

Close examination of the spikes caused by the McFarland compensation shows 
that wherever the spike occurs, the absolute value of the difference in the prediction of 
the aircraft state between two adjacent iterations is abnormally large (Fig. 4.1), much 
greater than the maximum of the corresponding absolute difference of the undelayed 
aircraft state. This fact leads to the definition of a criterion to decide if a spike will occur. 
If y and y p are the undelayed and predicted aircraft states, let 

md =Maxf =1 (|y[z] — y[z-l]|), 

where k indicates the current iteration, and let d p = y p [k] - y p [k - 1] , then, if 

r,=oj>=n (4.1) 

md 

a noticeable spike will occur (p varys with the time delay). For a delay up to 200 ms, 
p = 2.5 works well. If there is a noticeable spike, recalculate the prediction using 

*,[*] = *,[*-!] + % ( 4 - 2 ) 

and the spike will be reduced. Fig. 4.2 shows an example. 

Pilots flying the VMS at the NASA Langley Research Center commented that 
their perfonnance was much better with this algorithm. Nonetheless, it is worth 
mentioning that this simple spike -reduction algorithm is only an expedient measure. 
From Fig. 4.2, it’s obvious that the algorithm still causes prediction error, and more 
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PHI, rad 


importantly, the spike reduction algorithm reduces the phase lead in the spike- 
concentrated area. 


Delayed and Compensated Roll Angle with McFarland Filter 
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Fig. 4.1. Spikes caused by the McFarland compensation 
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McFarland Compensation of PHI with/without Spike-Reduction 
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Fig. 4.2. Spikes are reduced 

4.2. Frequency Domain Least Squares Method to Design McFarland Predictor 1 

In deriving the three coefficients of his predictor with sinusoidal tuning, 
McFarland used the system bandwidth, but he did not include a system dynamics model. 
For this reason, the McFarland compensator may work well for some systems with 
certain types of input, but it may not yield satisfactory results for other systems. The 
problem can be minimized if the three coefficients are designed by taking a system 
dynamics model into account. One way to achieve this is to make sure the frequency 


undelayed 

Compensation: Spike-reducted 


' If the coefficients are not determined with the sinusoidal tuning, the predictor is no longer technically a 
McFarland predictor, however, it is still referred to as a McFarland predictor for convenience in this 
discussion. 
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characteristics of the compensated system are as close to that of the undelayed system as 
possible, based on some defined criterion. The least squares method can be used for this 
purpose. As an example, a McFarland filter, designed with the frequency least squares 
method for the system depicted in Fig. 3.3, is presented. The design criterion is to 
minimize the cost function 

I = \ i \H d (Z, ) H McF (z,)- H 0 (z, )f (4.3) 

4- i=l 

where H d (z), H McF (z) and H 0 (z) are, respectively, the pulse transfer functions of the 
delayed system, the McFarland filter and the undelayed system, given by 

(z) = z| (*)][>„ (s)] [H td (5)]| (4-4) 


H McF (z) = 


1-z- 1 ^ 

V l + Z J 


+ 


{b 0 +b x z 1 +b 2 z 2 ) 




(4.5) 

(4.6) 


where z = e jT<> . The pulse transfer functions are used because the McFarland filter is 


discrete. In Eq. (4.4) and (4.6), the operation Z\ — [...]> is used to obtain the 


1-e 


-Ts 


l ^ J 

discrete transfer function using a zero-order hold (ZOH), with a sampling period T. 
H ac (s ) , H op (s) and H td (s) are given in Eq. (2.11), (2.10) and (2.4), respectively. 

Because the output of H ac ( 5 ) is the displacement (roll angle), and the McFarland filter 
uses the velocities, the differentiator s is pre-multiplied to change the output of H d [z ) 
to velocity. By choosing a suitable number N of frequency points 0 ) i within the pilot 
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working frequency range, the cost function (4.3) is minimized, and satisfactory results 
can be achieved. The final result is 



Although the cost function (4.3) seems to put constraints on the magnitude only, it 
actually minimizes the errors in both the real part and the imaginary part, so that the 
phase error is minimized also. The Bode diagram (Fig. 4.3) and step response (Fig. 4.4) 
of the compensated system with the McFarland filter designed with this method also 
shows that both the gain distortion and phase error are decreased compared with the 
ordinary McFarland filter. 
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Fig. 4.3. Frequency responses of compensated systems with McFarland filters 
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Fig. 4.4. Step responses of compensated systems with McFarland filters 
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4.3. Time Domain Least Squares Method to Design McFarland Predictor 

The frequency domain least squares method introduced in the previous section is 
suitable for designing a McFarland predictor in a simple simulation system similar to Fig. 

3.3. However, in a high fidelity, piloted flight simulator, H op ( 5 ) is replaced with a 

human pilot, and the flight dynamics are much more complicated than H ac (5) and may 

not be expressed as a simple transfer function. The frequency least squares method does 
not apply in this situation. Nevertheless, because the McFarland filter involves only the 
displacement and its velocity, which are available in the flight simulator, they can be used 
to design the McFarland filter coefficients that best fit the simulation condition according 
to similar quadratic criteria. A sample of roll axis displacement and velocity data 
recorded during an offset approach flown in the VMS at the NASA Langley Research 
Center, is plotted in Fig. 3.8, and is reused here to derive this predictor. Defining the roll 
angle and its velocity as y and v , respectively, the time domain quadratic cost function is 
given by 

1 = \Y,[_yd( k ) + K v ( k ) +b x v ( k - x ) +b 2 v ( k - 2 )- y( k Y\ ( 4 - 10 ) 

^ k = 1 


with k being the iteration index, and y d the roll angle y delayed by t d , which may be 
obtained from y using the Pade approximation. Minimization of the cost function I in 
the last equation results in a left pseudo-inverse as given in Eq. (4.7), with different 
matrices A and b . 


A r 


v(0) v (l) v(2) v(3) ■■■ v(N) 

0 v(0) v (l) v(2) ■■■ v(V-l) 

0 0 v(0) v ( l) ■■■ v(V-2) 


(4.11) 
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(4.12) 


j'W-j'rfW 

y(x)-yAx) 


Compensations with McFarland Predictor & LSF 3-V Predictor 



Fig. 4.5. Roll angle, velocity and roll stick of a simulation 

After the three coefficients b 0 , b t and b 2 are designed, calculate the predicted 

displacement with Eq. (3.11). A comparison of the McFarland compensators designed 
with the sinusoidal tuning and the least squares method, on the recorded data plotted in 
Fig. 3.8 is shown in Fig. 4.5. (In this case 192 ms prediction was applied.) With the least 
squares design, the spikes are significantly reduced, and the high frequency artifacts in 
the roll angle peak areas caused by the sinusoidally-tuned McFarland compensator are 
also removed. However, the calculation workload is substantially increased, and another 
serious drawback of the least squares methods is the three coefficients are not available 
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until the end of the simulation. Thus, this method is referred to as an offline least squares 
method. 


4.4. Adaptive Predictor 

The time domain least squares method discussed in the previous section can be 
implemented, for each iteration, so that the predictor coefficients are updated instantly. 
This involves calculating the pseudo-inversion (Eq. (4.7)) in each iteration, requiring 
intensive computation and storage. However, these calculations can be carried out in a 
recursive manner such that the current update uses only the current data and the results of 
the previous calculation. In this manner both redundant calculation and the storage of 
large quantaties of data are avoided. The coefficients are updated with the process so that 
the large prediction error caused by the constant filter coefficients is reduced. A 
schematic diagram of the adaptive predictor is illustrated in Fig. 4.6. 


Aircraft 

States 



Fig. 4.6. Structure of the adaptive predictor with the Kalman estimator 

To conveniently compare the adaptive predictor with the McFarland predictor, 
three consecutive steps of velocity will be used as an example. Therefore, using Fig. 4.6, 
the compensated aircraft state y c is given by y c = y d +b 0 v + b l v_ 1 + b 2 v_ 2 , where y d is the 
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delayed aircraft state, and v , v_ x and v_ 2 are the velocities of y in three consecutive 
iterations. Minimizing the quadratic loss function given by either Eq. (4.10) or 

l = \E[y-y c ] 2 (4.13) 

where operator “E” refers to the expectation. (Dividing Eq. (4.10) by N yields Eq. 
(4.13). The same results are obtained by setting the derivatives to zero. Eq. (4.13) was 
designed to conveniently introduce the ODE given in Eq. (4.30)). By making use of the 
Kalman matrix inversion theorem, the final recursive least squares method is given as 


P(*) = [l M -K(OfM] p (*-l) (4.14) 

K(*) = P(*-l)j(*)[l + J r (*)P(*-l)J r (*)]■' (4.15) 

e(*) = e(*-i)+K(*)[v(*)- (A-)-j r (*)e(*-i)] H-i6> 

where 

e(i-)=[MO MO M*ff (4.n> 

which gives the three coefficients of the adaptive predictor, and 

j r (*)=[v(*) v(k- 1) v(k- 2)] (4.18) 


is a vector consisting of the three consecutive velocities. The algorithm starts with 

P(M = (j(Mj r (M) 1 and 0(^o) = p (^ o )f {K){y( k o)-y d {K)) > where K 

corresponds to the first time when the quantity j r (k 0 ) j(£ 0 ) is nonsingular. Notice that 

since the quantity inside the brackets is a scalar, a matrix inversion is avoided, and the 
algorithm is considerably simplier than the original left pseudo inverse given by Eq. (4.7). 
This is the widely used Kalman filter algorithm. Because the coefficients are updated 
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each iteration, the prediction error is also reduced significantly. Simulations show that the 
new adaptive predictor substantially reduces the high frequency gain distortion and 
spikes caused by the McFarland filter. The three-step recursive least squares algorithm 
(Eq. (4. 14) - (4. 16)) can be formulated in a compact form as: 

§(*) = 8(*-l)+ , lW [y(k)-y d (k)-i T (k)~0{k-\)\ (4.19) 

ZiMfM 

i=h 

with ^j(z)j r (/) = P^(k) = P^(k-l) + j(k)j r (k). 

'=* 0 

The algorithm given by Eqs. (4.14)-(4. 16) allocates even weight to the data 
available no matter how old they are in the data history. A forgetting factor X(0<X<l) 

would allow the designer to assign larger weight to the more recent data so that they 
contribute more than the older data in the least squares algorithm. The new cost function 

1 k 

with the forgetting factor becomes / = — A k ~‘ 

2 t=k 0 

of the matrices P and K are changed to 

P(*) = [I M -K(*)j r (*)]p(*-1)/A (4.20) 

K(*) = P(*-l)j(*)[A+J r (*)P(*-l)r(*)]'' (4.21) 

The recursive fonnula for updating 0 is unchanged (the same as Eq. (4.16). An optimal 
forgetting factor range has been determined by trial-and-error and an algorithm with a 
forgetting factor in this range demonstrated a reduction of the predicting error. 

The recursive least squares method can be simplified for different approximation 
algorithms. The simplification comes from avoiding updating the matrix P , since 


( v (/)-y f . (/)) , and the update formulas 


73 



updating the matrix P dominates the computing effort for a large number of iterations. 
The first approximation considered was the Kaczmarz’s algorithm. But because the 
recursive least-square algorithm updates the current estimate 0 r (k) based on the previous 

estimate 0(£-l) and the new measurement y(k) - j 7 (k)(){k) , which contains 
information only in the direction of j r (k) in the parameter space, Kaczmarz proposed the 

nonnalized projection algorithm that minimizes 0 (/t)- 0 (/t) subject to the constraints, 
y(k) = j r (k)0(k) . The cost function for the Kaczmarz’s algorithm (the normalized 
projection) is / = |0(£)-0(£)|~+a[y(£)-j r (£)0(£)] , which may be considered as a 
function of the vector variable 0(/t) with a , the Lagrange multiplier, as a parameter. 
Taking derivatives with respect to 0(/fc) and a and invoking the stationary results in the 
Kaczmarz’s algorithm, also called the normalized projection algorithm, 0(/fc) becomes 

8(t)=e(t-i)+ . r( i ^ ) [j-(t)-^(0-f(t)e(t-i)] (4.22) 

Comparing the Kaczmarz’s algorithm with the original least squares algorithm given in 
Eq. (4.19) shows that the denominator of the Kaczmarz’s algorithm is changed from a 

k 

matrix ^ j(/) j 7 (/) to a scalar j r (k)j(k). This is where the simplification comes in. 

i=k 0 

Two other similar approximation least squares algorithms, with scalar denominators are 
the Stochastic Approximation algorithm and the Least Mean Square (LSM) algorithm 

e(0=e(fr-i)+-r-^ — [v(*)- (*)-j r (*)0(*-i)] (4.23) 

Zf (OKO 
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(4.24) 


0(k) = 0(k-l)+^[y{k)-y d (k)-j T {k)«{k-l)] 

The formula of the stochastic approximation (SA) algorithm ((4.23)) looks similar 
to the one-step basic least squares algorithm given in Eq. (4.19), with only the summation 
element being changed from a matrix j(/) j 7 (/) to a scalar j ( / ) j (/) . Inspired by the 

similarity, the stochastic approximation algorithm may also be rewritten in three steps: 
first Eq. (4.14) is rewritten as 

P(*) = [l-j(*)K(*)]p(*-l) (4.25) 

and the remaining two steps are exactly the same as Eqs. (4.15) and (4.16). If the 
exponential forgetting factor is applied to the three-step stochastic approximation 
algorithm, Eq. (4.20) becomes 

P(jfc) = [l-j(Jfc)K(Jfc)]P(Jfc-l)/A (4.26) 

and the formulas for the second and third steps are the same as Eqs. (4.21) and (4.16), 
respectively. Finally, the stochastic algorithm with a forgetting factor is formulated as 

0(*) = e(*-l)+— [y(k)-y d (k)-j T (k)0(k- 1)] (4.27) 

i=* o 

Offline tests of the five adaptive prediction algorithms, namely the Kalman filter 
algorithm, Kalman algorithm with a forgetting factor, the Kaczmarz’s algorithm, the 
stochastic approximation algorithm and the least mean square algorithm, on recorded roll 
angle data are illustrated in Fig. 4.7 and Fig. 4.8 (which is a blowup of Fig. 4.7) 

It is obvious from Fig. 4.8 that the stochastic approximation works best among 
these five adaptive algorithms. The next chapter will show the superiority of the 
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stochastic approximation by using three criteria. The stochastic approximation algorithm 
was chosen to implement the adaptive predictor in the final piloted tests. 


Comparison of Compensations with Four Adaptive Predictors 



Fig. 4.7. Adaptive compensations applied to the roll angle using different algorithms 
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Fig. 4.8. Zoom of Figure 4.7 

The recursive least squares methods have long been applied in system 
identification. In the past, the simplified recursive least squares algorithms, such as the 
stochastic approximation algorithm, were favored over the basic Kalman filter algorithm, 
because the former were more useful in real time when the processing speed of the digital 
computer was low. However, the simplified algorithms usually give biased identification, 
and they are vulnerable in the presence of noise. Therefore, as new generations of 
computers appeared with much higher processing speed, the basic algorithm became the 
preferred method once processing time was not a primary concern. Nevertheless, the 
stochastic approximation method is more suitable to application in a fight simulator when 
compensating for transport delay than the Kalman filter algorithm for three reasons. 

First, the application here is not system identification, rather it is a method to 
design a predictor aimed at providing accurate phase lead with small error. The input 
variables (the aircraft state and its velocity) contain little noise because they are results of 
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/'/ • 

• Ideal 
Kalman 

- 

• 

• 

- Kaczmarz 

— - - Stochastic 

- 

/'/ 

• 

Least mean square 

- 

l l # l 

1 1 1 




77 


the aerodynamic computation, rather than measurement, and one is simply the integration 
of the other. 

Second, the processing time is a concern. In order to minimize the transport delay, 
it is required that all the aircraft states in different coordinate systems be calculated in one 
frame after the operator control input is sampled. Unfortunately in the Langley VMS 
these computations consume almost the whole frame, leaving very limited time available 
for implementing the compensation. 

Third, the stochastic approximation algorithm yields the least predicting error 
(most accurate phase lead with least gain distortion) among all these recursive least 
squares algorithms. The explanation is given as follows: 

The basic requirement for a predictor of this type, i.e., using three consecutive 
iterations of velocity to predict, is b 0 +b l +b 2 = t d . This guarantees that the prediction is 

close to the time delay t d if the velocity change is not abrupt (low frequencies). But as 

analyzed in Chapter 2, if the velocity changes quickly (which frequently happens in real 
flight simulations), the filter causes a large error when the differences among the three 
coefficients are large. The benefit of using three steps of velocity instead of one is that a 
weighted average of three past values is less likely to cause a spike than just one past 
value. Another requirement of a three-velocity predictor is that the differences among the 
three coefficients must be small. A first choice might be b 0 -b l = b 2 = t d /3, but this is not 

necessary because the velocity changes irregularly. This is why an online recursive 
update of coefficients is employed. The coefficients may vary from simulation to 
simulation, but ideally they are as close to each other as possible. The smaller their 
difference, the better the prediction of the filter at high frequencies. Investigation shows 


78 



that, of the previously introduced recursive least squares methods, the stochastic 
approximation algorithm always generates the least difference among the coefficients. 
Table 4.1 shows the converged filter coefficients calculated with live different algorithms 
averaged across 16 sets of simulation data in the roll axis. Only the stochastic 
approximation method yields three coefficients that are close to each other. Note that the 
least mean squares algorithm alone gives all-positive coefficients in this case, but it has 
been found in other simulations that the coefficients have alternative sign changes. 

Table 4.1. Three coefficients calculated with different methods (7 rf =0.192s) 


Algorithm 

b 0 

k 


Sum 

McFarland 

Filter 

2.8613 

-5.2342 

2.5650 

0.192 

Kalman 

Filter 

-0.0030 

-0.5894 

0.7811 

0.189 

Stochastic 

Approximation 

0.0525 

0.0571 

0.0789 

0.188 

Kaczmarz 

Algorithm 

-0.2064 

0.0503 

0.3485 

0.192 

Least Mean 
Square 

0.0284 

0.0700 

0.1066 

0.205 


The algorithm expressed in Eq. (4.23) belongs to a large family called Stochastic 
Approximation Algorithms. Two other prominent versions are the Saridis and Stein’s 
algorithm and the Kwatny’s form . A general expression of the stochastic approximation 
algorithm is 

9(q=9(t-i)+ed(q[^(q-^(t)-j r (qe(t-i)] (4.28) 

where the scalar coefficient e k must be monotonically decreasing. Specifically Robins 
and Monro show that if 
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(4.29) 


e k > 0, Urn s k =0, V £ k = °° , and V £ k < °° 

k ^°° k~t ti 

the algorithm Eq. (4.28) converges. For the algorithm given in Eq. (4.23), 

1 

£ k =~ k 

Zf0')j(0 

i = k o 

where j r (A-) = [v(A-) v{k- 1) v(£-2)]. In a flight simulation, the velocity does not 

k 

change exponentially, so the coefficient £ k = ^ j r (z) j (/) meets the conditions given in 

i=k 0 

Eq. (4.29), and hence the algorithm Eq. (4.23) is indeed a stochastic approximation 
algorithm. Conversely, in the Kaczmarz’s algorithm £ k = j r ( / ) j (/) , and in the least mean 


square algorithm £ k is a constant, both of which do not meet the lim £ i = 0 requirement. 


Therefore, these are not stochastic approximation algorithms. 

The ODE that characterizes the asymptotic behavior of the stochastic 

23 

approximation algorithm given in Eq. (4.28) is 


= E(jf)-E(yf)0 


(4.30) 


where E is the mathematical expectation. The right side of the first equality of Eq. (4.30) 
is the negative gradient of the cost function, indicating that the stochastic approximation 
algorithm may be interpreted as a stochastic gradient descent algorithm. Though the 
gradient of the cost function in Eq. (4.13) with respect to 0 is unknown, the gradient at 

the current sample of -j[y(k)-y rf (&)]- j r (k)0(k)j is just 


j(k) y(k)-y d (k)- f (k)0(A:-l) , which is the dynamic term from Eq. (4.28). Because 
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the quantity inside the brackets y(k)-y d (k)-j T (k)Q(k-\) is a scalar, from the 
asymptotic ODE, 0 ( A: ) is in the direction of \ {k) . This means the modifying term of 
0(k + l) is in the direction of j ( k ) . And since the recursive algorithm starts with an 
initial zero vector, it is logical that 0 is in a direction close to the average of j . Because 
f(/) = [v(/) v(i- 1) v(z — 2)] , and the average velocity does not change much within 

two iterations, it follows that the three elements of 0 (or the coefficients of the adaptive 
predictor) do not differ much either. This demonstrates why the stochastic approximation 
algorithm gives the best prediction. 

4.5. A Practical State Space Compensator with a Reference Model 

As stated in section 3.3, although the Sobiski/Cardullo predictor shows some 
desirable advantages, it also has limitations in its implementation that prevent its practical 
application in a flight simulator. The basic state space prediction equation 
x p = Ox + TB« requires a linear time invariant (LTI) system, while the simulator flight 

dynamics are usually nonlinear, time-variant and coupled in different degrees of freedom, 
and are frequently not available. Instead of being expressed in state space equations, they 
are often expressed in coupled non-linear differential equations. However, employing an 
aircraft reference model in the predictor algorithm can solve this problem. This results in 
the development of a novel practical state space predictor, which is discussed as follows. 

4.5.1. Basic Implementation 

A reference model is an approximate linear aircraft dynamics model that is used 
to form the predictor states from the operator control input and the aircraft states, as well 
as to provide the state transition matrix and the convolution integral. This reference 
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model is used so that the state space prediction equation x p = Ox + *FB// may be applied 

in a flight simulator when the aerodynamics are not available. Fig. 4.9 illustrates this 
approach, where x a is the aircraft state vector (including the aircraft displacement, 

velocity and acceleration, etc), x is the fdter state vector and x p is the predicted fdter 
state vector. Because the aircraft state vector x a includes the aircraft state y to be 
predicted and its velocity and acceleration, the filter state vector x is also in terms of y . 
Therefore, the predicted filter state x p calculated with x p = Ox + TBz/ contains the 
predicted information of y . Then using the matrix C to retrieve y from x p , the desired 
prediction is achieved. 



Fig. 4.9. Structure of the state space compensator using a reference model 

Four 4 th -order reference models were tested. The first two models give the 
relationship between the pitch angle and the roll angle, respectively, and the 
corresponding stick deflections of a fixed wing jet flying at an altitude of 30,000ft and an 
airspeed of 430 kn ots. They will be called Model I and Model II (Eq. (2.11)). The other 
two models are for a large commercial transport, in the pitch mode, one for cruise, and 
one for landing. They will be called Model III and Model IV. These four models share 
the same general fonn of the transfer function 
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(4.31) 


b 2 s 2 + b v s + b 0 
s 4 + a 3 s 3 + a 2 s 2 + a x s + a 0 

The observable state space matrices of this general model are given in Eq. (4.32). 

Selection of the observable canonical form is made because the output is desired to be the 

first state variable. (In other words, y = x 1 ; and after the predicted filter state vector x is 

obtained, the predicted aircraft state is just the first element of x , i.e., y = x pl .) The 

expressions of the four state filter variables are directly derived fromJ , and 

[y = Cx + Dw 

the result is given in Eq. (4.33). Note that [xj x 2 x 3 x 4 ] r =x is the predictor state 


vector. This is an artificial state vector because, aside from the first element jq , the 


remaining elements do not exist in the simulator and have no physical significance. They 
are formed from the reference model. 


A = 


—a 3 

1 

0 

o' 


" 0 " 

-a 2 

0 

1 

0 

,B = 

b 2 






~Cl\ 

0 

0 

1 


_~ a o 

0 

0 

0_ 


_K_ 


, C = [1 0 0 0], D = 0 


(4.32) 


*i=y 

x 2 =y + a 3 y 

x 3 = y + a 3 y + a 2 y - b 2 u 

C t 

x 4 = 6 0 J 0 udt 


(4.33) 


Note that, in addition to Eq. (4.33), alternative formulas for calculating the state 
variables x, , x 3 and x 4 exist. Namely, x 2 = [ {-a 2 x x +x 3 +b 2 u]dt , 

JO 


x 3 = f {-a x y + x 4 +b 3 u\ dt and x 4 = x 3 + a l y-b l u . Though the final prediction, when 

J 0 

using the alternatives is only slightly different, these alternative equations are not used 
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because they require an extra integration, they include the jerk tenn x 3 , which introduces 
high frequency artifacts. 

In most aircraft simulations, the aerodynamic model is processed first, and then 
the accelerations and velocities of the vehicle are calculated in the body frame. These 
velocities and accelerations are then transformed to other necessary frames. For the visual 
display, the cueing channel for which the time delay compensation is designed, the 
topodetic coordinate frame is normally used (occasionally the geodetic frame , or north- 
east-up frame is used). In this coordinate frame, the six axes are roll, pitch, yaw, altitude, 
latitude and longitude. Prior to this study, the accelerations in these six axes were not 
normally available. The formulas to calculate these six accelerations are given below 
without their derivations (the derivations are in Appendix F). 


<i 


'1 

tan 0 sirup 

tan 0 coscp 

P 

0 

~ 

0 

coscp 

-sin <p 

q 

V 


0 

sin (p / cos 6 cos (p / cos 0 

r 


where (p , 0 and !// are the roll, pitch and yaw angles (Euler angles), and p , q and r are 


the three angular accelerations in the body frame. The other three accelerations are 


X 


u 


<i> 

/ 

= p 

A B2E 

V 

H 

+ 

0 

h 


w 




(4.35) 


where X , / and h are the latitude, longitude and altitude, and u , v and w are the three 
translational accelerations in the body frame, and the matrices P E2f , : and T are 


P r = 

1 B2E 


cos 6 cosy/ 
cos Osin (// 
sin 6 


sin cp sin 0 cos y/ - cos cp sin y/ cos cp sin 0 cos y/ + sin cp sin y/ 
sin (psin 6 sin yr + cos (/) cos (// cos tpsin 6 sin (// - sin (/) cos (// 
-sirup cos 0 - cos (p cos 6 


(4.36) 
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(4.37) 


«11 

oc n 

a 3l 

a 2\ 

0-22 

a 32 

a 3\ 

a 32 

a 33 


with the elements being given in Table 4.2. 

Another problem is on the input u in the six axes, since there are only four 
control inputs: the roll stick, the pitch stick, the rudder pedal and the throttle. For the roll 
and pitch compensation, the input is the roll stick and pitch stick, respectively. For yaw, 
both the roll stick and the rudder pedal contribute, but the input is chosen to be the rudder 
pedal because its control is more direct. Altitude is changed indirectly by the pitch angle 
or throttle, thus the pitch stick can be employed as its input. Longitude is affected by the 
heading, hence the rudder pedal may be used as its input. 


Table 4.2. Elements of the matrix T 


Element 

Expression 

Ot\ i 

v(cos <j) sin 6 cos i//+sin </>sin y/) + w(cos <p sin y/ - sin </>sin 6 cos y/) 

| OC\2 

-u sin 9 cos yr + v situ/) cos 9 cos y/+w cos (f> cos 9 cos yr 

a \3 

- u cos 9 sin yr -v [sin <p sin 9 sin yr + cos <j>cosy/) + w(sin </> cos y/ - cos </> sin 9 sin y/) 

l oc 2X 

v(co5 <j)sin9siny/-sin </>cos y/) - w(co5 <j)cos y/ + sin (j) sin 9 sin y/) 

\ 0-22 

—u sin 9 sin yr + v sin (j) cos 9 sinyr + wcos tp cos 9 sin y/ 

a 23 

u cos 9 cos y/+v ( sin </>sin 9 cos yr - cos <j) sin y/) + w(sin (p sin yr + cos <p sin 9 cos yr ) 

0-3 1 

-v cos (p cos 9+wsin<p cos 9 

oc i2 

u cos 9 + vsin<p sin 9 + wcos tpsin 9 

a 33 

0 


85 





















Comparison of State Space Predictor with Four Reference Models 



Fig. 4.10. Comparison of the state space predictors with four reference models 

After the filter state vector x is available from Eq. (4.33), apply the state space 
prediction formula x p = Ox + TBw to calculate the predicted filter state vector, and 

finally calculate the predicted aircraft state with y p = Cx p (D = 0). Compensated data, 

from the state space filters using the four formerly introduced reference models on 
recorded roll angle data are illustrated in Figs. 4.10 and 4.11 (an enlargement of Fig. 
4.10). They show that the state space filter with reference Model IV achieves the best 
prediction among the four reference models. The next chapter will present a detailed 
comparison of the effectiveness of the four reference models. Therefore, Model IV, or the 
landing model of a large commercial transport, in pitch has been chosen as the reference 
model for the state space compensation in the final tests. 
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Fig. 4.11. Zoom of Figure 4.10 

The reference models used in the state space compensation are not necessarily 4 th - 
order. Two 3 ld -order models 25 were also evaluated. The fonnulation of the filter state 
vector with the 3 ld -order reference model is similar to that with the 4 th -order, and hence 


only the result is listed here. The general transfer function and the state space expression 
of the model, and the filter state vector are, respectively 




b 2 s^ +b l s + b Q 
s' +a 2 s 2 + a { s + a 0 



~ a 2 

l 

o' 


b 2 

A = 

-a. 

0 

1 

, B = 



r a o 

0 
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K_ 


0], D = 0 


(4.38) 


(4.39) 
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(4.40) 


= y 

< x 2 =y + a 2 y-b 2 u 
x 3 = [ {-a 0 y + b 0 u\ dt 

JO 

There are also alternative expressions for some elements of the predictor state 
vector, but only the formulation that resulted in better compensation is shown. Although 
compensation achieved by the state space filter with a 3 ld -order reference model is not as 
good as with the 4 th -order models (the 3 rd -order reference model produces greater phase 
error but smaller gain distortion in the prediction), it is still better than the McFarland 
compensation. A detailed comparison among all these reference models will be presented 
in the next chapter. 

4.5.2. Simplification and Essence of the State Space Compensator 

Calculating the four predicted filter states as given by x = Ox + 'FB// involves 

many matrix operations. However, what is really needed is the predicted aircraft state y p 

given by y p = Cx , and because C = [l 0 0 0] , y is just the first element of x p . 

Therefore, calculation of the last three elements of x p is not necessary, and this shows 

that the algorithm can be simplified. The simplification is given in Eq. (4.41), where 

and Y] are elements of matrices <I> and *F (Eqs. (3.20) and (3.21)). By substituting the 


expressions of the four elements of the filter state vector in Eq. (4.33), the final simplified 
state space compensator is obtained as Eq. (4.42) 



y p =Cx p =(C<S>)x + {C'PB)u 
= [1 0 0 0]< 

— (j\ ]X[ + (j\ 2 X 2 "I" 03*L 014^4 Y \ U 

c T 

y p =(^ll+^2«3+^3«2)j ; + (<A2+^13«3)>’ + ^133 ) + (^l-tt)w+^oJ 0 udt ( 4 - 42 ) 

This is the essence of the state space filter! When compared with the previous 
compensator, this shows that, while the previous compensators use three consecutive 
steps of velocity to predict, the state space filter uses the current velocity, acceleration, 
the control input and its integral to predict. 

Likewise, the state space predictor based on the 3 rd -order reference model is 
simplified, and the final result is similar to Eq. (4.42): 

y P ={</\ i +^2«2)v+^ 2 j-^ 3 «oJ 0 ydt+iVi-fai^u+AA J 0 udt ( 4 - 43 ) 

Eq. (4.43) implies that the state space predictor may also be interpreted as a 
general PIP controller with two modification terms on the control input u. 

The prediction of the state space filter depends solely on the five coefficients in 
Eq. (4.42) or (4.43), which are functions of the time delay and the reference model. 
Therefore, an algorithm based on the state transition matrix and convolution integral with 
a reference model can be used as a design tool — to design the coefficients of the 
compensator. Because the reference model is time invariant, these coefficients are 
constants that may be calculated offline. Therefore, during each iteration of the 
simulation, only five multiplications and four additions are required — computation is 
simplified significantly. Because only the first rows of the matrices O and Y are needed 
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to calculate the predictor coefficients, calculating the remaining rows is unnecessary. 
This is significant in applications where a time-variant reference model is required, and 
the matrices <t> and *F must be updated every iteration. The state space filter may be 
implemented as depicted in Fig. 4.12 and Fig. 4.13, where the 90 ms of delay in the 
visual system is the baseline transport delay detennined in the VMS at the NASA 
Langley Research Center (artificial delay was also added for the experiments to 
determine the general applicability of the predictor to compensate other delays), and the 
coefficients are given in Table 4.3 and Table 4.4, respectively. 



Fig. 4.12. Simplified state space compensator using a 4 th -order reference model 

If the reference model is of 4 th -order, the coefficients of the last two terms of the 
control input u (c 3 and c 4 ) are so small compared to the first three terms without u that 
they may be neglected. In other words, a good reference model attenuates the 
contribution of the second tenn of \(t + t d ) =e A,J x(t) + ^ + r)dr, on which the 

state space predictor is based. The trivial contribution of the input u makes it easier to 
justify approximating the future input with the current input than Sobiski’s assumptions. 
On the other hand, the small contribution of the control input u is desirable because its 
high frequency jumps are attenuated. As the time delay gets longer, the coefficients c 3 
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and c 4 become larger so that more high frequency content is added to the prediction, 
resulting in a larger predicting error. Note also that the coefficient c 0 of tenn y is unity, 


a property the McFarland predictor also possesses ( y p ( k ) = y(k) + b 0 v + b l v_ l +b 2 v_ 2 ). 

Table 4.3. Coefficients of the state space predictor shown in Fig. 4.12 


Coefficient 

Expression 

For 4 th -order Model I, 
f,=186 ms 

c o 

01 1 "F 012^3 "F 013^2 

1 

c, 

02 "F 013^3 

0.1850 

C 2 

013 

0.0159 

c 3 

W\ — 013^2 

4.0052e-04 

c 4 

014^0 

3.2020e-05 
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Fig. 4.13. Simplified state space compensator using a 3 ,fI -order reference model 

Conversely, if the reference model is of 3 ld -order, the contribution to the 
prediction from the pilot control input u is greater from the integral tenn ( c 3 > c 2 ), and it 

may explain why the 3 ld -order reference model cannot achieve as good a prediction as the 
4 th -order reference model. The coefficient of y is no longer unity, and the longer the 
time delay, the farther it deviates from unity. This is a desirable property because as the 
time delay gets longer, the future y will resemble the current y less, and therefore the 
current y should contribute less to the prediction. 
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Table 4.4. Coefficients of the state space predictor shown in Fig. 4.13 


Coefficient 

Expression 

For 3 ld -order model A, 
t d = 186 ms 

c o 

011 T 012^2 

0.9949 

C 1 

012 

0.1750 

C 2 

1 

Is 

o 

-7.3938e-004 

C, 

ty\ 012^2 

0.0036 

C A 

013^0 

7.3938e-004 


4.5.3. State Space Predictor with a Discrete State Transition Matrix 

Chapter 3 introduced the discrete format of the state space predictor, and the 

discrete state transition matrix and convolution integral matrix (Eqs. (3.24)-(3.26)). A 

prediction algorithm based on the discrete state space filter using a reference aircraft 

model has also been tested in a manner similar to the continuous one. In changing the 

prediction algorithm from Eq. ((3.19) to Eq. (3.24), the fonnation of the filter state vector 

is also changed — it must be formed with the discrete state space equations: 

Jx(& + 1) = Gx(£) + Hm(£) 

| v(k) = Cx(£) + Dw(£) 

The general pulse transfer function corresponding to the continuous transfer 
function Eq. (4.31), its state space matrices and the discrete filter state vector are given as 
follows 


, b 3 z 3 +b 2 z 2 + b 1 z + b 0 

™ AC \ z ) 4 3 2 , , 

z + a 3 z +a 2 z +a l z + a 0 


(4.44) 


G = 


-a 3 

1 

0 

o' 


h 

-a 2 

0 

1 

0 

,H = 

b 2 

-a. 

0 

0 

1 


~ a o 

0 

0 

0) 


A_ 


C = [l 0 0 0], D = 0 


(4.45) 
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(4.46) 


x, (k) = y(k) 

x 2 (k) = -a 2 x l (k- l) + x 3 (k-\) + b 2 u (k-l) 

< 

x 3 (k) = -a.\X\ (k- l) + x 4 (k-l) + b x ii (k- 1) 
x 4 (k) = -a 0 x j (k-l) + 6 0 zz (k-l) 

Because the integrator in the continuous state space equation becomes the delay 

operator in the discrete equations, the derivative or integral terms are replaced with the 

past values of the aircraft state. Lack of aircraft velocity and acceleration information in 

prediction is the primary drawback of the discrete state space filter. Inclusion of the 

current and several past values may result in the derivative terms indirectly (due to 

numerical differentiation) introducing errors. Therefore, the discrete state space filter is 

expected to be inferior to a continuous one that makes use of the derivative terms directly. 

Application of the discrete state space predictor to the recorded simulation data has 

proven this. This is more obvious from the simplified algorithm. Corresponding to Eq. 

(4.42), the simplified algorithm with a 4 th -order reference model is given by 

y{k + l) = g n y(k)-(a 2 g l2 + a l g 13 + a 0 g l4 )y(k-l) + g n x 3 (k-l) 

+ g u x 4 (k - 1) + (b 2 g l2 + b lgl 3 + b 0 g l4 ) u (k - 1) + h t u ( k ) 

where g ij and /z, are the elements of the matrices G and H given in Eqs. (3.25) and 

(3.26). There are more contributions from the control input u , and conversely, the states 

x, (k-l) and x 4 (k- 1) have to be calculated recursively with the 3 rd and 4 th equations of 

Eq. (4.46), which results in error accumulation. Finally, as stated in Chapter 3, the time 
delay to be compensated with the discrete state space filter must be an integer multiple of 
the frame cycle. After considering all these factors, the discrete state space filter was not 
used during the final simulation tests. 
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4.5.4. Relationship Between Prediction and Reference Model 

Section 4.5.3 mentions that the four 4 th -order reference models give slightly better 
compensation than the 3 ld -order reference model, and the large commercial transport 
landing pitch model achieves the best compensation compared to the other three 4 th -order 
reference models. This raises some questions: which factors of the reference model make 
the most difference in compensation? Does the order of the reference model affect the 
compensation? What is the relationship between the prediction quality and the reference 
model? This section is a summary of the research designed to answer these questions. 


State Space Compensations with Two 3rd-order Reference Models 



t, s 

Fig. 4.14. State space compensations with two 3 rtl -order reference models 

The 3 ld -order reference model that achieved the best results is not an aerodynamic 
model, instead it was borrowed from a book on control theory (for convenience call it 3 ld - 
order Model A). Another 3 ld -order model (3 ld -order Model B) was also chosen, which 
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takes the same form of transfer function as given in Eq. (4.38). The state space 
compensations ( t d =0. 1 92s) using these two reference models on the same recorded roll 

data that was used for the other compensators are shown Fig. 4.14. 

The 3 rd -order model B gives a completely wrong prediction. Checks of the state 
transition matrix <t> and the convolution matrix do not reveal the problem. However, 
if the five coefficients c 0 — c 4 of the simplified filter shown in Fig. 4.13 are checked and 
compared with those of Model A, the problem becomes clear. These coefficients 
(C=0.192s) are listed in Table 4.5. The coefficients of the 4 th -order Model I, and the 
McFarland filter are also included for comparison. 


Table 4.5. Coefficients of different compensators for C=0.192s 


Filter 

Proportion 

Derivative 

Integral 

Input 

Integration 
of input 

c o 

C 


C 3 

c 4 

3 rd -order model A 

0.9946 

0.1803 

-0.0008 

0.0039 

0.0008 

3 rd -order model B 

0.4604 

0.0480 

-0.5299 

0.5659 

0.5299 

4 th -order model IV 

1 

0.1909 

0.0169 

4.3972e-004 

3.5149e-005 

McFarland 

1 

0.192 

(Z>, +Z>, +b 2 ) 





In Table 4.5, only the state space predictor with model B does not work. For other 
models, which work well, the coefficient of the proportional term is either unity or very 
close to unity, the coefficient of the derivative tenn is equal to or very close to the time 
delay, and all other coefficients are much smaller. This shows that the derivative 
extrapolation is either the only compensation or dominates the compensation. But this is 
not true for the 3 ld -order model B: its coefficient of the proportional tenn is far less than 
unity, and the derivative compensation contributes even less than the integral 
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compensation and the terms in the control input. This explains why the prediction quickly 
drifts away from the undelayed aircraft state. 

The values of the coefficients of the simplified state space compensator have 
direct relationships with the reference model. Table 4.4 shows that the coefficients are 
functions of the elements of the matrices O and Y , which are in tenns of a ' s , b's and 
t d . Deriving these expressions involves the evaluation of the matrices O and ¥ . The 
definition of the state transition matrix O is an infinite Taylor series 


O = e 


Atj 


I 

z=0 


(A t d ) 


il 


(4.48) 


Because this series is always convergent, it is usually approximated by truncating it to a 
finite series. The number of terms required to get satisfactory approximation depends on 
the matrix A . For the five reference models that work well, five terms seem to be 
sufficient. Therefore, 


<p = e A '" - 1 + A t d + | ( At d f + ^ ( A t d ) 3 + ^ ( At d ) 4 


(4.49) 


In fact, the last term contributes little to the final result. Another way to compute 
the state transition matrix O is the exact method making use of the Cayley-Hamilton 
theorem 26 . That algorithm is given by 

n - 1 

<J> = e Atj = Y J a A l ( 4 - 5 °) 

;=o 

where n is the order of the matrix A , and a 0 - a n _ x are the solutions of the following 
linear equations 

\ = e 1 '*' 1 . U = l-») (4.51) 

i = 0 
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with the Aj ’s being the eigenvalues of the matrix A , provided that all the eigenvalues are 


distinct. The analytical expressions of a 0 - a n _ x with/without repeated eigenvalues are 
given in Appendix A. For the five working reference models, it has been found that 


, 1 2 1 3 

OCq 1, cc x t d ,cx 2 t d , oc 2 t d 
2 6 


(4.52) 


Note that only the 4 th -order models have a, . The approximations in Eq. (4.52) 


can be verified by the Taylor expansion of e Xjt, ‘ : e Xj,d =1 + (^/ (/ ) + — + (^/ </ ) /*!+••• For 


t d up to 0.3 seconds, and the absolute values of the real part of the maximal eigenvalue 
around 1, the fifth term [Ajt d ) /24 is only 3.3750e-004, and the higher order tenns are 
even smaller. Therefore, the scalar Taylor series can be truncated to four terms, or 

e Xjtd + {A J t d) + [A j t d ) 12 + [Ajt d ) 16. 

Then Eq. (4.51) becomes 

1 + [Ajt d j + (Ajt d j” / 2 + {Ajt d ^ / 6 ~ (Xq + cc x Aj + cc 2 A~ + (4.53) 

Comparing the coefficients of both sides gives Eq. (4.52). In this sense, the state space 
compensation may be viewed as a Taylor series extrapolation in the state space form. 

The convolution integral matrix approximated by making use of Eq. (4.49) is 
given in the last equation of Eq. (4.54). Substituting the elements ^ and i// ] of the 

matrices <I> and T calculated with Eq. (4.49) and Eq. (4.54) into the expressions of the 
five coefficients listed in Table 4.4 (for state space predictor based on a 3 -order 
reference model; application on a 4 th -order reference model will be discussed later), the 
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final expressions of c 0 ,c x ,...,c 4 are given in Table 4.6. Due to the approximations in Eqs. 
(4.49) and (4.54), these expressions are also only approximate . 

Y = p' e A(,d ~ T) dz 

Jo 

= (a>-l)A _1 (4.54) 

= t d + — A + — A 2 + — A 3 
2 6 24 

Table 4.6. Approximate expressions of the coefficients of the simplified state space 
predictor with a 3 rd -order reference model 


Coefficient 

Expression 

c o 

2 , 2 

i 2 a^ci 2 a 0 ^ ^ 0^2 ci^a 2 ~i 4 

2 ti+ 6 ‘ 1+ 24 h 

c l 

, d 2 2 ^2 ^1 ^3 . 2a { a 2 ci 2 +a 0 4 

d 2 d 6 d 24 ! ‘ 


2 

*0 f 2 . ^ 1^2 A . ^ 0^2 * 0*1 A 

Irf • Cj 1 Tj 

2 6 24 

C 3 

/ 1 2 A / 1 \ 

1 2 Cl] t. Cl 9 ^14. 1 2 4 r 

-t d +— t d +— — -r d b x + -t] — U d b 0 
{2 d 6 d 24 d J 1 [6 d 24 d J 0 

C 4 

/ 2 \ 

1 2 d] 2 a*-) ci 1 4 i 

— L ^+— ~t d b 0 

[2 d 6 d 24 d J° 


Further examination of the expression of c 0 , indicates that the third and fourth 
tenns are much smaller than the second term for a time delay shorter than 0.3s. And the 
second term, a x must be less than 0.5 for a delay of 0.2s or less than 0.25 for a delay of 
0.3s for c 0 to be unity within 1%. For example, for model A, a ] =0.3050 results in 
c 0 =0.9946 which would be satisfactory; for model B, a x = 62.1314 results in c 0 =0.4406, 
which is much less than unity. With a similar analysis, a 2 must be less than one to make 
c, fairly close to t d . Thus, the values of a x and a 2 are critical to the applicability of a 3 rd - 


98 

















order model as a reference model. The restrictions on a ] and a 2 indicate that there must 
also be some restrictions on the eigenvalues of the reference model because, for a 3 ld - 
order model, the following relationships between the three eigenvalues \ ^ and a x , a 2 
are held 

f A, +/!.,+ A-, = ~~ct 2 

4 - ^ 2 (4.55) 

[W + + W ~ a \ 

For a stable model, all eigenvalues must have negative real parts. Then Eq. (4.55) implies 
that \A t \ <a 2 <1 and A I A / < a, < 0.5 . Therefore, the absolute values of all eigenvalues must 

be at least less than unity to make the reference model work. Model A meets this 
requirement, but model B does not. The eigenvalues of model B are located much farther 
from the imaginary axis of the s-plane than model A, and therefore, model B responds 
much faster than model A (Fig. 4.15). In other words, the bandwidth of model B is much 
larger than that of model A. 

Therefore, it appears that a good reference model can be formed by merely 
choosing suitable a's and b's in Table 4.6 so that the five coefficients have desirable 
values. Although it appears that the reference model does not have to match the flight 
dynamics, in actuality this is not true. In fact, the reason the 3 ld -order model B does not 
work for the roll angle recorded from the VMS is because it has totally different dynamic 
properties than those of the aircraft model running in the VMS. In other words, the reason 
that Model B does not work as well as Model A does is that Model A has frequency 
characteristics close to the aircraft model used in the simulation from which the data was 
collected, but the Model B does not. Likewise, the four 4 th -order reference models and 
the 3 ld -order model A work well because their bandwidths are close to the bandwidth of 
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the large commercial transport model running on the VMS. This can be verified by 
simply comparing the bandwidths of the 3 ld -order models A and B with those of the four 
4 th -order aerodynamic models, as shown in Table 4.7. The bandwidth of the 4 th -order 
model II is not available because it has a pure integrator (Type I system, Eq. (2.1 1)). 


Table 4.7. Bandwidths of six reference models 


Model 

3 rd -order 

4 th -order 

A 

B 

I 

II 

III 

IV 

Bandwidth 

(rad/s) 

0.6069 

8.0561 

0.6523 

NA 

0.1036 

0.0324 


Step Responses of Two 3rd-Order Reference Models 



t, s 


Fig. 4.15. Step responses of two 3 rd -order reference models 

The bandwidth of the 3 ld -order model B is much too high compared to the other 
aircraft models. One way to determine whether the model bandwidth plays an important 
role in compensation quality would be to reduce the bandwidth of model B gradually, and 
apply it to the state space prediction, and see how the compensation changes. The results 
are demonstrated in Fig. 4.16. Investigation shows that satisfactory compensation can be 
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achieved when the bandwidth of the reference model is below 1 rad/s. This proves that 
the bandwidth of the reference model must be close to that of the simulated vehicle 
dynamics. 

Originally, it was assumed that the damping ratio of the reference model made a 
significant difference in compensation. Some of the reference models that do not work 
well are over-damped, whereas the aerodynamics of the large commercial transport are 
under-damped. A study of the impact of varying the damping ratio was conducted, 
similar to the study in which the bandwidth was varied, and that study shows that the 
damping ratio of the reference model does not make a significant difference in the state 
space compensation. 


State Space Compensation Using Model B with Varying Bandwidth 



Fig. 4.16. State space compensation using model B with varying bandwidth 

rd 

With a derivation similar to that of the 3 -order reference model, the five 
coefficients of the simplified state space predictor with a 4 th -order reference model as 


101 


depicted in Fig. 4.12 were obtained in tenns of the parameters of the model transfer 
function. These coefficients are listed in Table 4.8. By comparing Table 4.8 with Table 
4.6, it becomes obvious that there is no term in t] . Because t d is smaller than t] for 
? rf <ls, the proportional coefficient c 0 , when using a 4 th -order model, is closer to unity for 
the same value of a x , and the derivative coefficient c ] is closer to t d for the same value 
of a 2 . Furthermore, because the coefficients a x and a 2 are higher order functions of the 
eigenvalues when using a 4 th -order model (compare Eqs. (4.55) with (4.56)), even though 
the eigenvalues of a 4 th -order model have magnitudes close to those of a 3 rd -order model, 
and the a t and a 2 of the 4 th -order model are smaller (the absolute values of the aircraft 
model eigenvalues are usually less than 1), this makes c 0 and c, closer to unity and t d , 

respectively. Thus, the 4 th -order reference model is expected to achieve better prediction 
than the 3 rd -order reference model, and the comparison between compensations with the 
3 rd -order model A and 4 th -order model IV proves this. 

(A l A 2 + A l A 3 + A l A , 4 + A 2 A 3 +A 2 A 4 +A 3 A, 4 =a 2 (4 56) 

yA^A^A^ + A 2 A 2 A , 4 + A 4 A 2 A 4 + A 4 A 2 A 4 = — 
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Table 4.8. Approximate expressions of the coefficients of the simplified state space 
predictor with a 4 th -order reference model 


Coefficient 

Expression 

c 0 

1 f 3 | a \ a i ~ a 0 ^4 

6 d 24 d 


^ a 2 {2 | ~ a \ ^4 

d 6 d 24 d 

C 2 

1 2 

1 9 ci'j q ci-> a, a 

-t] — 3 ~r d + — -~t d 

2 6 24 

C 3 

( a, 3 a, 4 V b n 4 
—t] — 3 -r d b 2 +—t d 
U 24 d ) 2 24 d 

C 4 

f I n \ 

, l d l d u o 

\6 24 J 


The comparison between these two reference models may not be convincing 
because the 3 rd -order model A is not an aircraft model. So, a 3 rd -order model was formed 
by reducing the 4 th -order model, such that the reduced-order model shares similar 
frequency characteristics with the original model. Then this reduced-order model was 
used to carry out the state space compensation. The compensation error was indeed 
considerably greater than that of the original 4 th -order model. 

It has been shown that the 4 th -order reference model is superior to a 3 rd -order one. 
But, what about a 2 nd -order model, or a model of order higher than 4? If the reference 
model is of 2 nd -order, the simplified state space filter is given as 

y P = oc 0 y + a x y + ( 1 - a 0 ) u (4.57) 

n — 1 

where a 0 and a x are solutions of the coupled equations i,2), 

(=0 
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where A 0 and \ are the two eigenvalues of the model. Eq. (4.57) is inferior to a reduced 

state space predictor with a 4 th -order reference model because less system information is 
used for prediction. On the other hand, if a 5 th -order reference model is employed, the 
filter state vector will contain either a high frequency jerk component (derivative of the 
acceleration) or triple integration of the aircraft state, which is likely to introduce 
computational artifacts. In short, the 4 th -order reference model is the best choice. 

As stated previously, the bandwidth of the reference model plays a major role in 
the compensation quality of the state space predictor. The bandwidth of a model is 
affected by its poles and zeros, but not by its gain. The gain of the reference model also 
influences the compensation quality. From Tables 4.6 and 4.8, the coefficient c 3 of the 
control input u is a linear function of the model gain. A large gain of the reference model 
will amplify the high frequency components in the control input u , and distort the 
prediction. However, the sensitivity of a reference model’s compensation quality to the 
variation of its gain depends on the individual model; some models are more sensitive 
than others. The effects of the gain of a reference model are much less significant than 
the effects of its bandwidth. 

The large commercial transport landing model in pitch (4 th -order model IV) has 
been found to work successfully on data recorded from the large commercial transport 
simulation runining in the LaRC VMS. If the vehicle dynamics were quite different, e.g., 
a fighter rather than a transport, one possible method of finding a suitable reference 
model with which to apply the state space predictor would be to use system identification 
techniques to get a model that closely matches the aircraft dynamics. During the course 
of a simulation run, the dynamics may change dramatically. For instance, the dynamics 
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before and after an offset maneuver may be quite different. In such cases, fuzzy logic 
may be employed to identify the mode. Unfortunately, applying these algorithms 
increases the computation workload greatly. These issues merit further investigation. 

In summary, the following points are conclusions of the study of the state space 
predictor, and the relationships between the compensation quality and the reference 
model: 

• The state space predictor can be simplified to a general PID filter; 

• The state space compensation may be viewed as a Taylor series extrapolation in 
the state space fonn; 

• The bandwidth of the reference model must be close to that of the simulated 
vehicle; 

• The damping ratio of the reference model does not make a significant difference 
in the state space compensation; 

• The 4 th -order reference model is the best choice; 

• The effects of the gain of a reference model are much less significant than the 
effects of its bandwidth. 
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5. Results of Theoretical Analysis 


In Chapter 4, two novel predictors for compensating the transport delay in a flight 
simulator are compared from a theoretical perspective: the adaptive predictor based on 
the Kalman filter algorithm, and the state space predictor using an aircraft dynamic 
reference model. Chapter 4 demonstrated that these two new types of compensators show 
improved performance over the McFarland predictor. It also mentioned that the stochastic 
approximation adaptive predictor achieves the best compensation among the five 
different Kalman filter algorithms, and the large commercial transport landing model in 
pitch works best as a reference model for the state space predictor. All these conclusions 
are made without sufficient verification while focusing only on the principles of the novel 
compensators. In order to give a convincing comparison among all those compensators, 
quantitative approaches are necessary to analyze the respective compensation results. 
Because different predictors of the same type (i.e., the five adaptive predictors, and the 
state space predictors using different reference models) may demonstrate slight 
differences in compensation quality, the quantitative methods must be sensitive enough to 
highlight the differences and give accurate results. 

This chapter starts by defining two metrics, which can be used to evaluate the 
compensation errors caused by different predictors. It then applies these metrics to the 
compensation results from offline tests, and finally quantitatively compares the 
compensation qualities among different predictors in terms of the two metrics. In the 
offline tests, compensations were applied to the aircraft state data recorded during a 
previously completed simulation study, generating predicted aircraft states. Offline tests 


106 



were implemented on a personal computer running MATLAB rather than the Visual 
Motion Simulator. No visual images were generated and no pilot was involved. The 
predicted aircraft states were used only for analysis. The last section investigates the 
sensitivity of the prediction errors with respect to the amount of time delay. 


5.1. Error Metrics 

Compensation error may appear as either phase error or magnitude error (gain 
distortion) or both. Fig 5.1 illustrates an example of phase error and magnitude error of a 
compensation applied to an ideal sinusoidal signal. 


Illustration of Phase Error and Gain Distortion of Compensation 



0 0.5 1 1.5 2 2.5 3 



0 0.5 1 1.5 2 2.5 3 



Fig. 5.1. Illustration of phase error and gain distortion of compensation 

The function of a predictor is to generate a phase lead designed to be equal to the 
transport delay in the succeeding subsystem so that the delay can be compensated by the 
phase lead. In offline tests, the transport delay may be simulated by simply delaying the 
state information by the desired amount. In Fig. 5.1, the “compensation” signal is the 
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result of the delayed prediction. The ideal compensation produces no error (either phase 
error or distortion) at all, as depicted in the upper subplot of Fig. 5.1, where the 
compensated signal lines up exactly on top of the undelayed one. In other words, at every 
discrete point, the values of the two signals are exactly the same. Define an error function 
as the sum of square difference between the compensation and the undelayed signal 

N r , 

£ = S[v c (z)-v(0] C 5 - 1 ) 

;= l 

where y and y c are the undelayed and compensated signals. This error function would 
be zero for the ideal compensation. 

Assume that the predictor brings no gain distortion, but the phase lead it provides 
is not exactly equal to the phase lag resulting from the time delay. This is shown in the 
middle subplot of Fig. 5.1, in which the compensated signal does not match the 
undelayed signal. The phase error causes the predicted value to not match the undelayed 
signal at most time points (except those intersections D and E). For example, when t=1.5s, 
the phase error is AC , and the corresponding magnitude difference is AB . Because there 
is no gain distortion, the greater the phase error, the larger the magnitude discrepancy if 
the phase error is less than k (A phase lag of n results in an opposite phase, representing 
the maximal magnitude error E. As the phase changes from 0 to n , E increases, but as 
the phase lag exceeds n , E starts to decrease). Therefore, if there is no gain distortion, 
either the phase error or the error function defined by Eq. (5.1) may be used as a metric to 
evaluate the compensation error. 

A real compensator usually results in both phase error and gain distortion, as 
shown in the dashed dot curve in the bottom subplot in Fig. 5.1. Evaluating the error 
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index E, as given by Eq. (5.1), between the dashed curve and the dashed dot curve, both 
phase and gain errors are present, but, the error index does not indicate which factor 
contributes more. For the transport delay compensation, the phase error usually is given 
more attention, and therefore it is necessary to define two metrics of error so that the 
effects of the two factors can be separated. 

One way to separate the two is to fit the compensated signal with a smooth curve 
(such as the thick dashed curve in Fig. 5.2). Determine the phase difference between the 
fitted curve and the undelayed curve, and define this phase difference as the phase error 
of compensation. Then calculate the error index E between the actual signal (dashed) and 
the fitted one (dashed dot) using Eq. (5.1). In this example, the undelayed signal is a 
known function. Because the fitted signal is also analytical, it is easy to determine the 
phase difference between these two analytical functions by simply drawing a horizontal 
line intersecting both curves, and the distance between the two intersection points is the 
absolute value of the phase error. 
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Roll Angle, Prediction and Polynomial & Sinusoidal Approximation 




Fig. 5.2. Roll angle, prediction and polynomial & sinusoids approximation 

In most situations, the real signal to be compensated is not an analytical function. 
In such situations, the undelayed signal also needs to be fitted with an analytical function. 
Because the actual aircraft state in a flight simulator is much more complicated than a 
sinusoidal signal, the fitting function is a combination of many simple components (a 4 th - 
order polynomial and 13 sinusoids of various magnitudes, frequencies and initial phases). 

5.2. Comparison of Predictors Based on Offline Tests 

Based on the analyses of the compensator used in the offline tests, in terms of the 
two error metrics defined in the previous section and the magnitudes of spikes on the 
compensations, this section will present a comparison of compensation qualities among: 

1) The McFarland predictor and four adaptive predictors. 

2) Five state space predictors using five different reference models (four 4 th -order 
models and one 3 rd -order model). 
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3) The McFarland predictor, the adaptive predictor using the stochastic 
approximation algorithm, and the state space predictor using the large commercial 
transport landing model in pitch. 

5.2.1. Comparison of the McFarland Predictor and the Adaptive Predictors 

These five predictors were applied to the aircraft roll angle data recorded from 16 
straight-in approach tests. The predictors were designed to compensate for 48, 96, 192 
and 288 ms (all integer multiples of the update period 16 ms) of transport delay, and in 
each time delay category, the amount of prediction (in ms) was averaged across the 16 
test runs. Table 5.1 gives the average predictions and standard deviations for each case. 

Among the four adaptive predictors, the least mean square algorithm 
demonstrates larger phase prediction error and a larger standard deviation than the other 
three. The basic Kalman filter, Kaczmarz and stochastic approximation algorithms show 
insignificant differences in phase prediction, whereas the stochastic approximation 
algorithm tends to have smaller standard deviation than the previous two when the 
transport delay is long (equal or greater than 192 ms). All four adaptive algorithms 
generate larger phase prediction error than the McFarland compensator, but apart from 
the least mean square algorithm, the differences are not significant. The standard 
deviations of the stochastic approximation algorithm and the McFarland predictor are 
very close, while the mean phase prediction error of the fonner is about 5 ms larger than 
that of the latter. 
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Table 5.1. Mean values & STD of the predictions with five 3-velocity predictors 


t d (ms) 

Mean Prediction (LEFT) and Standard Deviation (RIGHT) (ms) 

McFarland 

Predictor 

Basic Kalman 
Predictor 

Kaczmarz 

Algorithm 

Stochastic 

Approximation 

Least Mean 
Square 

48 

48 

1 

43 

1 

43 

1 

43 

1 

41 

5 

96 

96 

2 

93 

3 

93 

2 

93 

2 

90 

8 

192 

190 

5 

185 

8 

191 

15 

187 

6 

179 

18 

288 

283 

9 

271 

14 

275 

14 

278 

10 

266 

24 


Table 5.2 lists the compensation gain errors defined in Eq. (5.1) of these five 
predictors averaged across the same 16 test runs. Among the four adaptive algorithms, 
the basic Kalman filter algorithm causes the greatest gain error, and the stochastic 
approximation algorithm causes the least. Apart from the stochastic approximation 
algorithm, all the other adaptive algorithms show larger gain distortion than the 
McFarland predictor. As mentioned in Chapter 4, the McFarland predictor generates 
noticeably larger spikes than all the adaptive predictors, and this indicates the error 
measure in Eq. (5.1) is a more useful metric of gain error than the magnitude of spikes. 
The former is defined for a whole test, whereas the spikes only exist locally. In other 
words, greater total gain error does not mean larger spikes. If the transport delay is short, 
the gain error of compensation by the stochastic approximation algorithm is even greater 
than the McFarland predictor, but as the time delay becomes longer, the former displays 
smaller gain distortion. 
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Table 5.2. Gain error index of the McFarland and four adaptive predictors 



Gain Error Index 

t d (ms) 

McFarland 

Basic 

Kalman 

Predictor 

Kaczmarz 

Stochastic 

Least Mean 


Predictor 

Algorithm 

Approximation 

Square 

48 

0.0416 

0.3818 

0.3731 

0.2705 

0.4763 

96 

0.7735 

2.5586 

2.3050 

1.0363 

1.9108 

192 

17.6686 

37.2153 

33.1492 

13.9038 

22.8639 

288 

74.5940 

171.2643 

152.5877 

59.4277 

95.3306 


The motive for developing the adaptive predictors is to reduce the large spikes 
and gain distortion of the McFarland predictor. The analyses given above of Table 5.1 
and Table 5.2 demonstrate that the stochastic approximation algorithm is superior to the 
other three adaptive algorithms. Therefore, an adaptive predictor using the stochastic 
algorithm is the best among the five in this group. 

5.2.2. Comparison of Five State Space Predictors 

In Chapter 4, live aircraft reference models for the state space predictor were 
introduced, i.e., four 4 th -order models: the pitch model (Model I) and the roll model 
(Model II) of a fixed wing jet flying at an altitude of 30,000ft and an airspeed of 430 
knots, the cruise model of the large commercial transport in pitch (Model III) and the 
landing model of the large commercial transport in pitch (model IV), and a 3 rd -order 
model. The gain errors, phase predictions and standard deviations of compensation by 
these live predictors for the four transport delay cases are listed in Table 5.3 and Table 
5.4, respectively. While Model III introduces the least gain distortion, the difference in 
terms of gain error between Models III and IV (the two large commercial transport 
models) are negligible; the two fixed-wing jet models introduce very similar gain error 
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(the pitch model brings slightly smaller gain error than the roll model), but the gain 
distortion is obviously larger than that of the two large commercial transport models; and 
the 3 ld -order model introduces the largest gain distortion. 

Table 5.3. Gain error index of the state space predictors with five reference models 


t d (ms) 

Gain Error Index 

4 th -order 
Model I 

4 th -order 
Model II 

4 th -order 
Model III 

4 th -order 
Model IV 

3rd-order 
Model A 

48 

0.0071 

0.0080 

0.0062 

0.0062 

0.0542 

96 

0.1378 

0.1573 

0.1063 

0.1079 

0.7374 

192 

4.2312 

4.5971 

2.8558 

3.0095 

9.9674 

288 

28.5543 

30.4164 

16.1325 

17.5617 

44.4455 


Table 5.4. Mean values & STD of state space prediction with five reference models 


t d (ms) 

Mean Prediction (LEFT) and Standard Deviation (RIGHT) (ms) 

4 th -order 
Model I 

4 th -c 

Moc 

rder 

el II 

4 th -order 
Model III 

4 th -order 
Model IV 

3rd-c 

Mod 

urder 
el A 

48 

48 

1 

48 

1 

48 

1 

48 

1 

48 

1 

96 

94 

2 

95 

2 

95 

2 

96 

2 

95 

3 

192 

175 

5 

184 

10 


4 


4 

186 

8 

288 

236 

10 

266 

31 

275 

8 

287 

5 

273 

16 


Model IV yields the least mean phase prediction error and the least standard 
deviation (from 16 test runs). For delay up to 192 ms, the mean phase is zero, and for a 
delay of 288 ms, the phase error is only 1 ms. Models II, III and the 3 ld -order model 
introduce very similar mean phase errors (about 1 frame for a delay of 288 ms), though 
the standard deviation of Model II is significantly greater than the other two when time 
delay is 288 ms; Model I introduces the greatest mean phase error among these five 
models. 
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Therefore, in both gain distortion and phase error, the large commercial transport 
landing model in pitch (Model IV) works best compared with the other four reference 
models. In general, the 3 rd -order model is inferior to the 4 th -order models. 

5.2.3. Comparison of the McFarland, the Adaptive and the State Space Predictors 

From the comparisons in the previous two sections, the stochastic approximation 
algorithm is the best adaptive predictor, and the large commercial transport landing 
model in pitch is the best candidate as the reference model for the state space predictor 
among the five models. Thus, these two novel predictors were chosen, to be implemented 
along with the McFarland compensator in the final piloted tests, in order to compare them 
in a human-in-the-loop simulation. In this section, the perfonnance of these three 
compensators will be compared with offline tests. Though the data are available in Tables 
5.1 - 5.4, the mean phase prediction and gain error specifically for these three predictors 
were extracted and are listed in Table 5.5. As an example, Table 5.3 illustrates the phase 
compensations of these three. 
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Table 5.5. Mean values & STD of the predictions of three predictors 


t d (ms) 

Mean Prediction (LEFT, ms) and Gain Error (RIGHT) 

McFarland Predictor 

Adaptive Predictor 

State Space Predictor 

48 

48 

0.0416 

43 

0.2705 

48 

0.0062 

96 

96 

0.7736 

93 

1.0363 

96 

0.1079 

192 

190 

17.6686 

188 

13.9038 

192 

3.0085 

288 

283 

74.5940 

278 

59.4277 

287 

17.56175 
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Fig. 5.3. Phase lead generated by three types of predictors (^=192 ms) 

The state space predictor tends to introduce significantly less gain distortion than 
the other two. It causes less phase prediction error than the adaptive predictor, and 
achieves slightly better phase compensation than the McFarland predictor for long time 
delays. It seems to be superior to the adaptive predictor, but it introduces moderately 
larger spikes than the latter. The phase prediction error of the adaptive predictor is larger 
than that of the McFarland predictor, but not by much. When the delay is longer than 192 
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ms, the adaptive predictor significantly reduces the high frequency artifacts, especially 
the spikes generated by the McFarland compensator, resulting in the reduction of gain 
error. Although the McFarland predictor shows the least phase error, which is determined 
using a smooth fitting curve, the actual phase error in areas of large gain distortion may 
be even larger than that of the other two predictors. 

5.3. Sensitivity Analysis 

As the transport delay to be compensated gets longer, there is no doubt that both 
the gain distortion and phase error increase. But how fast the errors increase for different 
compensators, or in other words, how sensitive they are to changes in time delay, is still 
an interesting problem. This section is an investigation of the sensitivity of three 
compensators: the McFarland predictor, the adaptive predictor and the state space 
predictor. 

In order to complete this analysis, the group of compensators was implemented, 
and tested on the 16 sets of aircraft roll angle data recorded from piloted simulations. For 
this test, the transport delay was varied from 48 ms to 288 ms, and the frame time was 16 
ms. The phase prediction and gain error were calculated, and the results are summarized 
in Table 5.6, which is actually an expansion of Table 5.5. 

The phase prediction of the state space predictor is insensitive to changes in time 
delay, whereas the McFarland predictor and the adaptive predictor are sensitive to 
changes in time delay to a similar degree (from 0 to 5 ms for the former, and from 5 to 1 0 
ms for the latter, as the delay increases from 48 to 288 ms). In terms of absolute increase 
of the gain error, the state space predictor is the least sensitive to the time delay, and the 
McFarland is the most sensitive (the absolute gain error increases of these three 
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predictors are, in descending order, 74.4524, 59.1572 and 17.5556). But in terms of 
relative gain error increase, the adaptive predictor is least sensitive to the change of time 
delay, while the state space predictor is most sensitive. The relative gain errors increase 
by 1791, 220 and 2833 times respectively for the McFarland predictor, the adaptive 
predictor and the state space predictor. But, the most important metric is the absolute gain 
error increase, and the state space predictor is the least sensitive to the time delay, 
whereas the McFarland predictor is the most sensitive in both the phase and gain errors of 
prediction. 

Table 5.6. Mean predictions and gain error index of the McFarland predictor, an 
adaptive predictor and a state space predictor ( t d e [48, 288] ms) 


t d (ms) 

Mean Prediction (LEFT, ms) and Gain Error (RIGHT) 

McFarland Predictor 

Adaptive Predictor 

State Space Predictor 

48 

48 

0.0416 

43 

0.2705 

48 

0.0062 

64 

64 

0.1127 

93 

0.2786 

64 

0.0099 

80 

80 

0.3202 

77 

0.5211 

80 

0.0402 

96 

96 

0.7736 

93 

1.0363 

96 

0.1079 

112 

112 

1.6210 

109 

1.8812 

112 

0.2332 

128 

128 

3.0330 

125 

3.1277 

128 

0.4477 

144 

144 

5.1968 

141 

4.8614 

144 

0.7962 

160 

160 

8.2767 

156 

7.1684 

160 

1.3126 

176 

175 

12.4131 

171 

10.1465 

176 

2.0359 

192 

191 

17.6686 

187 

13.9038 

192 

3.0095 

208 

206 

24.1493 

202 

18.5569 


4.2808 

224 

222 

30.1740 

217 

20.4700 

224 

8.6872 

240 

237 

40.8535 

232 

31.0103 

240 

7.9571 

256 

253 

50.9008 

248 

39.0578 

255 

10.5185 

272 

268 

62.1538 

262 

48.4900 

272 

13.6921 

288 

283 

74.4940 

278 

59.4277 

287 

17.5618 
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In addition, as the time delay becomes longer, all three predictors tend to be more 
sensitive to the time delay. For instance, for the McFarland predictor, the gain error index 
increases only 0.071 1 as the delay increases 1 frame from 48 ms, but it increases 12.3402 
as the delay increases 1 frame from 272 ms. In other words, the prediction error does not 
increase linearly with the transport delay. 
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6. Conclusions and Future Research 


6.1. Conclusions 

Theoretical analyses in both time domain and frequency domain show that the 
effects of the transport delay depend on the simulation dynamics. Specifically, the higher 
the system bandwidth, the larger the loss of the phase margin, or stability. Meanwhile, the 
system suffers with the same amount of transport delay, and therefore, the system has less 
ability to tolerate the time delay. In a flight simulator system, higher bandwidth means 
the system is more vulnerable to pilot-induced oscillation and instability. 

Measurements were made with a device called SIMES to measure the transport 
delay in the Visual Motion Simulator (VMS) at the NASA Langley Research Center. The 
transport delays in the visual system, the motion system and the instrument system were 
measured to be 58, 56 and 40 ms, respectively, with the frame cycle being 16.7 ms. The 
average total transport delays from the pilot control input to the ends of these three sub- 
systems were determined to be 89.7, 87.7 and 71.7 ms, with the frame length of the 
simulation computer being 16 ms. Therefore, the average baseline transport delay of the 
visual cue is 90ms (rounded from 89.7ms) . This forms the basis for the time delay 
compensator designs for the final piloted simulation tests. 

Three basic criteria (phase error, gain distortion and computation complexity) 
were adopted to assess and compare the three prominent, previously developed time 
delay compensators — the lead/lag filter, McFarland predictor and the Sobiski/Cardullo 
predictor. This was done using theoretical simulations with the same aerodynamic model 
and pilot model as used for the analysis of time delay in Chapter 2. For the lead/lag filter, 
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the formula proposed by Crane to design the fdter pole was revised, resulting in better 
compensation. By applying the McFarland compensator to the recorded simulation Euler 
angles, large spikes in the compensated output resulting from its gain distortion have 
been found and detennined to be the primary cause of the poor compensation observed in 
the preliminary piloted tests. Investigation shows that the excessive gain distortion comes 
from the alternative sign changes of the three coefficients of the McFarland predictor, and 
from the constant coefficients, which do not change with the simulation process. For the 
Sobiski/Cardullo predictor, a discrete version was formulated which yields the same 
compensation results as the fonner in theoretical simulation. Some limitations in 
Sobiski’s implementation were found. 

A spike -reduction algorithm was proposed as an expedient solution to the 
annoying large spikes in the McFarland compensation. Both time and frequency domain 
least squares methods were formulated to design the coefficients of a three-velocity 
predictor. A novel adaptive predictor was developed using the Kalman filter algorithm 
that updates the predictor coefficients during the course of a simulation. Five versions of 
the adaptive algorithms were tested and compared, the stochastic approximation 
algorithm was chosen for the final piloted tests because it produces the smallest 
compensation errors. A stochastic approximation algorithm with a forgetting factor was 
also formulated. The analysis in Chapter 4 mathematically proves the viability of using 
the asymptotic ODE, and demonstrates why the stochastic approximation algorithm 
achieves the best compensation among all the five adaptive predictors. 

The first practical state space predictor to compensate the transport delay was 
developed by forming a group of predictor states and making use of a linear aerodynamic 
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reference model. This makes it possible to apply a state space compensator to a flight 
simulator in which the aircraft dynamics are nonlinear and time-variant. The state space 
predictor was significantly simplified from the original Sobiski/Cardullo implementation, 
and reduced to a general PID filter. Several 3 rd and 4 th -order dynamic reference models 
were tested and compared, and the large commercial transport landing model in pitch was 
chosen for the final piloted tests. A mathematical explanation of why this model is the 
best reference model for use in the state space compensation was discussed. While 
exploring the relationships between the state space compensation effectiveness and the 
reference model, the following important points were found: 

1) The state space compensation may be viewed as a Taylor series extrapolation in 
the state space form; 

2) The bandwidth of the reference model must be close to that of the simulation 
dynamics; 

3) The damping ratio of the reference model does not make a significant difference 
in the state space compensation; 

4) The 4 th -order reference model is the best choice in this case; 

5) The effects of the gain of a reference model are much less significant than its 
bandwidth. 

Theoretical analyses were conducted by applying the McFarland predictor and the 
two novel predictors to recorded Euler angles from past simulation tests, and evaluating 
their phase and gain errors. The adaptive predictor with the stochastic approximation 
algorithm produces slightly higher phase error than the McFarland predictor, which 
reveals about the same phase error as the state space predictor with the large commercial 
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transport landing dynamics as the reference model. The difference in phase errors of the 
three predictors is not significant, and therefore the gain errors play a more important role 
in deciding their relative effectiveness. The adaptive predictor yields a greater gain error 
than the McFarland predictor for the zero and 48 ms added delays (short delays), but 
shows a smaller error for the 96 and 192 ms added delays (long delays). The advantages 
of the adaptive predictor become more obvious for longer time delay. Conversely, the 
state space predictor results in substantially smaller gain error than the other two 
predictors for all the four delay cases. 

6.2. Suggested Future Research 

It would be worthwhile to study a frequency domain method to measure the 
simulator transport delay and phase lag, and compare the results with those obtained with 
the time domain method. Although some frequency domain data were collected with the 
SIMES while it generated sweep signals as control inputs to the VMS, they did not result 
in satisfactory results because aliasing was included in the aircraft EOM outputs. The 
frequency of the sweep signal was increased too quickly, and the lower frequency input 
did not last long enough. A high percentage of high frequency input signals caused the 
aliasing. This could be avoided by choosing a lower rate of increasing the frequency of 
the sweep signal. 

The effectiveness of the three-velocity adaptive predictor for compensating a 
short time delay may be improved by updating the coefficients in a more “intelligent” 
manner, which is more adaptable to the velocity changes, especially sudden and abrupt 
changes, such as the initiation of a maneuver. An exponential forgetting factor can only 
weigh the past data continuously and gradually, and therefore cannot solve this problem. 
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It would be worthwhile to develop some “intelligent” algorithms that can detect 
significant changes in the flight dynamics. Fuzzy logic and neural network techniques 
may be helpful to achieve the detection of these changes. Maneuver phase detection 
would be of great significance, and may be useful in applications other than transport 
delay compensation. 

The compensation quality of a state space predictor depends on the choice of the 
reference aircraft model. As is stated in Chapter 4 (4.5.4), system identification may be a 
good way to obtain an optimal reference model. To avoid adding extra delay resulting 
from a highly complex aircraft model, techniques to minimize the computational burden 
should be investigated. Because different phases of a maneuver may involve varying- 
order dynamics, system identification may also require maneuver phase detection. 
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Appendix A. Calculation of State Transition Matrix and Its Convolution 
Integral 


A.l. State Transition Matrix 

Computation of the state transition matrix d> and the convolution integral matrix 
*F of the state matrix A plays an important role in implementing the Sobiski/Cardullo 
filter and the aerodynamic reference model state space predictor. These matrices can be 
calculated directly or indirectly. The direct method makes use of the definition of the 
state transition matrix, and the indirect method makes use of the eigenvalues of the 
matrix A , or the inverse Laplace transform. 


A. 1.1. Direct Method 

By definition, Eq. (4.49), the state transition matrix is an infinite Taylor series 


O = e Xtj = Yj 


(AC)' 


(A.l) 


Z! 


Because this series is always convergent, the infinite series can be approximated with a 
finite number of terms. The number of tenns required depends on the criterion set for the 
accuracy. If the power of the matrix A is equal to or higher than its order n , this power 
can be calculated in tenns of those lower-order powers by making use of the Cayley- 
Hamilton theorem 


^(A) = X«,A'=0 (A. 2) 

/=o 

where a 0 - a„ are the coefficients of the characteristic equation of A . This shows that the 
A" can be calculated in terms of A, A 2 , .... A" 1 , i.e., 
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A. 1.2. Indirect Method 

The Cayley-Hamilton theorem, Eq. (A. 2), says that every square matrix satisfies 
its own characteristic equation. This leads to the indirect calculation of the state transition 
matrix in that the infinite series in Eq. (A.l) can be rewritten as a product of 0(A) and a 

polynomial in A plus a remainder, that is 
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0 = e Ai ' =P(A)0(A) + Z?(A) 


(A. 7) 


where the remainder R (A) is of lower order than 0(A) . The counterpart of Eq. (A.7) in 
tenns of a scalar variable is given by 

e x = P(x) 0(x) + R(x) (A. 8) 

From the Cayley-Hamilton theorem, 0( A) = 0 , then 


<I> = e 


Xtj 


R( A) 


: 0CqI + a x A + ... + oc n _ x A 


2—1 


(A. 9) 


The coefficients a lt ,a can be calculated by noting that every eigenvalue of A 


satisfies 0(A ) = 0 , substituting this into Eq. (A. 8) gives 


e ^ =R(A t ) = a 0 + a l A t + ... + oc n _ x A," 1 (A. 1 0) 

If the matrix A has all identical eigenvalues, the n equations in (A. 10) can be solved for 

the coefficients a n ,a a n _ x . Solving the coupled equation group involves inverting a 

Vandermonte matrix 


1 A , A , 2 ■■■ A ,'- 1 

1 A A 2 A '" 1 

1 A A - AA 


(A. 11) 


The solution is 


1 

if 


i 

i 

if 

II 

< 

A 

e Xl ‘ d 



1 

^ • 
1 


If there are repeated eigenvalues, with a multiplicity of m ( m<n 
more equations are needed, and they are given by 


(A. 12) 
), then (m- 1) 


127 



TJT^) = ^K + ^ 


,j = — 1 


(A. 13) 


If a pair of eigenvalues of complex conjugates exists, the resulting state transition 
matrix is real because the two complex equations can be replaced with two real equations 
obtained from the real and imaginary equalities. For example, for A = p± jq , the two 

complex equations e hd = R(A) = a 0 + a 1 A + ... + a n _ l A”~ 1 can be separated into two real 


equations 


e p ' d cos[qt d ) = real^OC 0 + a x (p + jq) + ... + a n _ 1 (p + jq)" 'J 
e p ' d sin(qt d ) — imaglcx 0 +a l (p + jq) + ... + a n _ 1 (p + jq)" 'J 


(A. 14) 


Another indirect approach for calculating the state transition matrix is to use the 


inverse Laplace transform,, which is given by 


' = e A ‘ d = X 1 {[^1 - A]' 1 } 


(A. 15) 


Inversion of [,s l - A] can be computed using the expansion of its adjoint as follows 


[si- A] -1 = 


adj(sl- A) 

Ui-aI 


(A. 16) 


where si - A is the characteristic polynomial: 


si -A =s"+a )I _ 1 s" ' +a n _ 2 s n z + ... + a { s + a 0 


(A. 17) 


The coefficients a 0 -a„ can be calculated with Eq. (A.4). The adjoint matrix may be given 


adj(sl- A) = Is" 1 +T„_ 2 s" 1 +... + T 1 S 1 +T 0 


(A. 18) 


where the matrices T, - T„ are given by Eq. (A. 5) 
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A.2. Calculation of the State Transition Matrix Integral 


By definition, the integral of the state transition matrix is given by 

y = (A. 19) 

If the matrix A is non-singular, it can be calculated by direct integration 

*F = A^ 1 [e Xtd -i) = A -1 (<D -I) (A. 20) 

But if A is singular, this method cannot be used because its inverse matrix does not exist. 
Matrix inversion can be avoided by substituting the definition of <D into Eq. (A.20), 
which then becomes 



(A. 21) 


The repetition of completing the factorial may be avoided by using a recursion algorithm 


Y = fJl + A^ l + X T A I + A%Tl + ... + fl + A 


(A. 22) 


Using the recursion in a manner similar to the implementation of Eq. (A. 6). The 
calculation of the matrix T is very similar to calculation of <t> . The power functions of 
the matrix A used when forming 0> with the direct method may be saved for later 


calculation of T . 


If the order of A is not high, the method given by Eq. (A.20) is simpler than the 
recursive algorithm (A. 22) when the matrix 0> is already available. However, if A is of 
high order, algorithm (A.20) becomes complicated, because inversion of a high order 
matrix is time consuming. In this case, use of the recursive algorithm is recommended. 
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A.3. Calculation of the Discrete State Transition Matrix 

The calculation of the discrete state transition matrix <1> , = G is straightforward, 
however, it can also be obtained using the inverse Z-transfonn, as follows 

0>, = G' = Z' 1 {[zl - G] 1 z] (A. 23) 

l J k—l 
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Appendix B. Complements to Sobiski/Cardullo Filter 

Sobiski and Cardullo first introduced the idea of using the state transition matrix 
to compensate the time delay in flight simulation. The state space predictive filter named 
after them bears many unique features for designing and implementing state space 
compensation. This section presents a brief description of the Sobiski/Cardullo predictive 
filter and some additional work done by the author of this report. 

B.l. State Space Compensation: Output Feedback and State Feedback 

Figure 2.3 is a block diagram of the flight simulation system, which includes the 
Sobiski/Cardullo compensator as it was implemented for this study. The transfer 
functions of the operator model, the aircraft dynamics and the Pade approximation of the 
time delay were first transformed to state space equations and output equations, which 
were then concatenated together, resulting in a single-input and single-output (SISO) 
linear time-invariant system, represented by 

f x = Ax + Bw 

(D = 0). 

[y = Cx 

Because the Pade approximation of the time delay is included, the output y has 
been delayed, and the state space prediction is to be employed to compensate for the 
delay. The matrices A , d> and T all contain information about the time delay, and are of 
10 th -order. The state space compensation may be implemented as depicted in Fig. B.l. 
The compensated system is described by the equations 

f x = Ax + B u 

(B.l) 

[v = CO>x + C , FBw 

The matrix D - C'FB , a scalar now (a one-by-one matrix), is not necessarily zero. 
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Fig. B.l. State space compensation in an open loop system 

The unit output feedback closed-loop system with state compensation is shown in 
Fig. B.2, and the state equation and output equation of the closed loop compensated 
system are 


Jx = (A-BGCO>)x + BGv 
|v = CO>Gx + C’FBGv 


(B.2) 


where G = (l + C*FB) ' is, as Sobiski and Cardullo defined, the feed forward gain of the 


closed loop compensated system. 
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Fig. B.2. Output feedback closed loop system with state space compensation 

The state feedback can also be used in a state transition matrix-based 
compensation. Actually the Sobiski/Cardullo predictive filter uses the state feedback, as 
illustrated in Fig. B.3. If the state feedback is used, the feed forward gain matrix becomes 

G = (l + K'FB) ' . The feedback matrix K is designed using a pole placement technique 
so that the closed loop poles of the compensated system consist of the closed loop poles 
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of the undelayed system and the real parts of the Pade approximation of the time delay. 
The first step is to obtain K fdhk using ordinary pole placement so that the homogeneous 

equation x = (A-BK^)x has the desired poles. The second step is to calculate the 

feedback gain K as follows: comparing x = (A-BK /(M )x with the first equation of 

(B.2) gives K fdbk = GKO , and the equations G = (1 + K'FB) and K fdhk - GKO 
together give 

K=K H „(®-«PBK / „y (B.3) 



Fig. B.3. State feedback closed loop system with state space compensation 

Because the real parts of the 2 nd -order Pade approximation are identical, the 
desired closed loop poles for the pole placement have repeated elements, and this 
prevents the MATLAB “place” function from working. A new MATLAB function has 
been written to achieve pole placement with repeated poles (Available in Appendix C.14, 
NASA CR— 2007-215096) 

The final Sobiski/Cardullo filter is a little different from Fig. B.3, and is actually 
an equivalent system obtained using the feed forward gain G . The equivalent system is 
shown in Fig. B.4. It can be proved that the systems in Fig. B.3 and Fig. B.4 share the 
same closed loop transfer function 
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(B.4) 


T c (s)_ K<D(sI- A) 1 BG 
R{s) 1 + K® (xl - A) -1 BG 



Fig. B.4. The Sobiski/Cardullo filter 

One advantage of the equivalent system shown is that the compensated system 
may be easily simplified to canonical fonn because the matrix *F in the initial state 
feedback system (Fig. B.4) is changed to the feed forward scalar gain G . However, these 
two systems are not exactly equivalent although the feedback y c is the same, the system 
state x is different, and hence the system output y is also different. Comparing the two 
systems shows that they are identical only if the feed forward gain G is unity. Because G 
is not unity when using the pole placement method, the system states are changed. In an 
aircraft simulation, the aircraft states may be used by different cueing systems. With the 
equivalent system, although the operator is fed back with the correct, compensated visual 
cue, the instrument cue and motion cue (if any) are not correct. This is a limitation of 
Sobiski/Cardullo compensation. 

When the system is of high order, e.g., above the 10 th -order, pole placement may 
sometimes cause problems. The state feedback matrix K can be designed using other 
approaches. Two frequency domain methods have been used: the five-point method and 
the least squares method. Both methods assume the operator working frequency does not 
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exceed an upper limit co 0 . In the five-point approach, five frequency points are chosen 

evenly from the pass band [0&> 0 ], and the frequency characteristic function of the 

compensated system is forced to equal that of the undelayed system at these points. It is 
given by 

K [yfy,I - A] ' B + *FbJ = H 0 (y'ry ( ) , (/ = 1 5) (B.5) 

where H 0 is the Laplace transfer function of the undelayed system. The five real part 

equalities and five imaginary part equalities can be used to solve the 10 elements of the 
matrix K . 

In the least squares approach, more than five frequency points are chosen, and the 
sum of the Euclidean norm of the difference between the frequency characteristic 
functions of the compensated system and the undelayed system at these points is 
minimized. In mathematical tenns, minimize the equation 

S = X K A] 1 B + 'Pb}-// 0 (y'ft) ) " (B.6) 

Z=1 

Then, setting the derivatives of S with respect to each element of the matrix K to zero 
results in 10 coupled linear algebraic equations that can be used to solve K . Analysis 
shows that the state feedback matrix designed with the least squares methods can achieve 
better compensation than with the five-point method and much better than with the pole 
placement method. Figs. B.5 and B.6 show the Bode diagrams and step responses of the 
state space compensation with the feedback matrix K designed by the least squares 
fitting (LSF) method for delays of 200, 400 and 800 ms. The compensation is perfect and 
the responses of the compensated system line up exactly on top of the uncompensated 
response. 
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Bode Diagrams of a System with LSF-Determined Gain S/C Compensations 



Fig. B.5. Bode diagrams with state space compensation with K by LFS 


Step Responses of a System with LSF-Determined Gain S/C Compensations 



Fig. B.6. Step responses with state space compensation with K by LFS 
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Another advantage of the five-point method and the LSF method over the pole 
placement method is that the gain matrix K , designed with them, makes the forward gain 

G = (l + K*FB) ' very close to unity, and because of this, the equivalent system is almost 

identical to the original system, without a noticeable change of system states. In the live- 
point method and the LSF method, the frequency characteristic function of the 
compensated system is forced to be close to that of the undelayed system. Hence the 
matrix D of the open loop compensated system ( D - C*FB ) is also forced to be close to 
that of the undelayed system, which is zero. Therefore, D - K'FB ~ 0 , and 

G = (l + K*FB) 1 ~ 1 . This explains why the feed forward gain is unity from the two 
frequency design methods. This is significant in the simplification of the compensation. 

B.2. State Observer for the State Space Compensation 

If the state space filter is used in another control system to compensate for time 
delay, and some of the states to be fed back are not measurable, then a state observer can 
be designed to make available those immeasurable states. Fig. B.7 illustrates a state space 
compensation system making use of a full order state observer. The state equations and 
the output equations of the control system and the predictor with the observer are, 
respectively 

x = Ax + B u 

<v = Cx (B.7) 

u = v-y c 

ji = A.,x + B«, + Lv 
[y = CO>x + CTO< 

where A c = A - LC . Define e = x - x then it follows from Eqs. (B.7) and (B.8) that 
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e = x-x 


(B-9) 


= A c e 



If all eigenvalues of A c are selected to be negative, Eq. (B.9) is an asymptotically 
stable error equation. The eigenvalues of A c control the speed of convergence. 

Experience indicates that a good design usually results if the continuous-time observer 
poles are selected to be a little farther to the left in the s-place than the desired closed 
loop state feed back poles (Brogan 1974). This technique is used here. The fact that the 
control system is now completely observable makes it possible to choose any eigenvalues 
for A c as desired. By eliminating the variable u (u = G(v-C3>x) ), the equations that 
characterize the entire system are obtained as 
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X 

= 

A 

-BGCO 

X 

+ 

BG 

.ij 


LC 

A e -BGCO> 



BG 


v = [ca>-c , FBG(ca>) o] 


+ C*FBGv 


(B.10) 


where G = (I + C*FB) ' . The observed states can be fed back with the gain matrix K 


designed with the two frequency methods introduced above in place of the matrix C . The 


poles of the observer and the poles of the feed back controller can be chosen 
independently because the poles of Eq. (B. 10) can be separated as follows 


A t (A) 


1,,/l-A BGKO> 

-LC I n A - A c + BGKO 


1,,/L-A + BGKO) BGKO> 

o M-A, 

|I„/l-A + BGK®||I„/l-A I 


=a c (1)a'(1) 


(B.l 1) 


As previously stated, using the feedback matrix K , which was designed with the 
two frequency methods, makes the feed forward gain G be close to unity and the 
equivalent compensation system identical to the original compensation system. This also 
makes it possible to adopt the equivalent system to the observer system, as illustrated in 
Fig. B.8. The state equation of the observer is changed from Eq. (B.8) to 

x = Ax + BGw + LC(x-x) (B.12) 


where the control law is also changed to 

= v- v = v-K3>x (B.13) 

From these, the governing differential equation of the error for the new observer 
becomes 


e = x-x = (A-LC)e + B(l-G)w 


(B.14) 
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Since G ~ 1, this equation reduces to Eq. (B.9). Therefore, the design of the observer 
gain matrix L may be undertaken the same way as for the observer in Fig. B.7. 

If some state variables can be measured, then a minimal order state observer may 
be designed to implement the state space compensation. 
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Appendix C. State Space Compensation in LTI Systems 


In the state space compensation, introduced in Appendix B, the Pade 
approximation of the time delay is added first so that the matrices A , and T all 
contain the information about the time delay. The filter generates the compensated states 
from the delayed states. In this manner the prediction is not visible. In this section, the 
prediction y is produced before the time delay is added, and then it is delayed to get the 

compensated output to use as feedback. This approach is illustrated in Fig. C.l, where the 
matrices A a , <t> and T , do not carry the information about the time delay. The transfer 
function of the 2 nd -order Pade approximation of the time delay is given in Eq. (2.4), and 
the A, , B, , C, and D, are given by 


A, = 


6 





1 


12' 



B,= 

t d 

12 

0 




0 

. c 2 





C,=[l 0], D, = 1 


(C.l) 


In whichever form is chosen for the state space equations, the matrix D is always 

unity because the Pade approximation transfer function is normal and its gain is unity. 
Then the delayed output can be calculated using 


1 

•>< 

1 


i 

> 

& 

O 

I 

1 

M 

& 

i 


B„ 

1 

j*!- 

i 


B,C fl A, 

I 

i 

+ 

0 


y d = K c, 


x, 


(C.2) 


The prediction is given by 


\K =A a X a +R a U 

] v = C 0> x +CPB u 

p a a a a a a 


(C.3) 
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The predicted state is output to the cueing system to compensate for the time 
delay. The compensated output may be calculated in a manner similar to that used for Eq. 
(C.2) as follows 


K 


A a 0 

X a _ 


B, 



_B,C a O A, 

_ X <_ 

+ 

bc.tb, 


y c = K*a c, 


\ a 

X, 


+ CJPB u 


(C.4) 


Finally, the output feedback closed loop state equation and output equation are 
given by Eq. (C.5), where G - (l + C„*F fl B u ) ’ is the feed forward gain. The 
compensation results shown in Figs. 4.14 and 4.15 were obtained with this algorithm. 


x „ 

X, 


A a B a C,0> a -B,C, 

BA(I-WJ3>„ A, -B ( C u 'F u B l C, 

X, 


B a 

B / C„T a B a 


y c = G[ c„o a c, 


(C.5) 


+ GC u 'P„B,v 



Fig. C.l. State space predictor in a linear time-variant system 

The matrix C a can be replaced with a matrix K designed with similar frequency 
methods described in Appendix B. With the gain matrix K , the feed forward gain 
becomes G = (l + K*F fl B a ) _1 . 

As shown in Figs. 3.13 and 3.14, this compensation based on prediction cannot 
compete with the Sobiski/Cardullo filter when compensating for long time delay. For 
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example, the compensation errors, which come from compensation based on prediction, 
are considerably larger then the errors arising from the Sobiski/Cardullo approach as the 
time delay approaches 800ms. In addition, it has been demonstrated that the longer the 
delay is, the farther the feed forward gain G deviates from unity. Remember that, as the 
time delay gets larger, the greater compensation errors mean the compensated system 
diverges from the undelayed system, and hence the matrix D c = CTB of the 

compensated system moves away from D fl of the undelayed system, which is zero. This 
also means that G — (l + C a 'F a B fl ) 1 deviates further from unity. Nonetheless, the 

prediction is significant. The state space predictor introduced in this section is the first 
step in the evolution from the Sobiski/Cardullo filter to the aircraft reference model state 
space filter applicable to a flight simulator. 
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Appendix D. Discrete State Space Filter and Time-Variant State Space 
Filters 


D.l. Discrete State Space Filter 

The discrete state space control system of the corresponding continuous control 

fx = Ax + Bu {x(k + \) = Gx(k) + Uu(k) 

system < may be expressed as < , where G = e 

[ v = Cx \y(k) = Cx(k) 

and H = (| 0 e dv |b The discrete state space predictor comes from the recursive 

formula for calculating the succeeding states in terms of the previously available ones 

x(£ + 2) = Gx(A: + l) + Hw(£ + l) = G 2 x(A:) + GHw(A:) + Htt(A: + l) 

x(k + 3) = G\(k + 2) + Uu(k + 2) = G 3 x(k) + G 2 Uu(k) + GUu(k + l) + Uu(k + 2) 


(D.l) 


k - 1 

x(k + l) = Gx(k + l-l) + Uu(k + l-l) = G'x(k) + G 2 Uu(k) + Y J G k ~ J ~'Jlu(j) 

j = 0 

By making the same assumptions on the input u as for the continuous state space filter, 
the future inputs can be approximated by the current input u(k), and the discrete state 

space predictor can be obtained from the last equation of (D. 1), as given by 

x(k + l) = <b d x( y k) + x P d u(k) (D.2) 


/-i 

where <& d (/) = G ; and x V d (I) = Z<M»h are the discrete state transition matrix and 

j = o 

its integration matrix. The predicted output is 

y(k+l) = C® d x(k) + C'¥ d u(k) (D.3) 

The unit output feedback discrete state space compensation may be implemented 
as shown in Fig. D. 1 . 
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Fig. D.l. Discrete state space compensation 


D.2. Continuous Time-variant State Space Filter 

The time varying continuous state space filter can be derived from the solution of 
the time varying continuous state equation x(t) = A(t)x(t) + B(t)w(t), which is given 
by 

x(?) = O(?,? 0 )x(/ 0 )+ f ®(?,r)B(r)w(r)dr (D.4) 

J / 0 

where ®(?,r) is the general state transition matrix defined as 


®(r,r) = I. 
d<J>(t, t) 


= A(t)&(t,r) 


(D.5) 


By making some assumptions about the control input u and the system matrix A , the 
continuous time-variant state space predictor is given by 

x(t+t d )= e A ^' d x(t)+ ^ e A ^ td ~ T) dr B [t)u{t] (D-6) 
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D.3. Discrete Time-variant State Space Filter 

The time-variant discrete state space filter comes from the solution of the time- 
variant state equation x{k + 1) = G(k)x(k) + H(k)w(k) , and is given as 


k - 1 
j=K 


x ( k ) = 0> ( k , k 0 ) x ( k 0 ) + £ <D ( k, j + 1 ) H ( j ) 1 1 ( j ) 


(D.7) 


with <i>(k,k 0 ) being the time -variant state transition matrix, which can be calculated by 


®(*oA) = I - 

<t>(k,k 0 ) = G{k- l)G(k-2)---G(k 0 ) 


(D-8) 


By making some assumptions about the control input u and the system matrix G , the 
discrete time -variant state space predictor is given by 

r i - 1 i 


x[k+l)= G [k)‘ x(k) + 


XG(k) 7 H(k) 


>(*) 


(D.9) 


j = o 


The state transition matrix G [k)‘ 
to be calculated for each iteration. 


/-i 

and its numerical integral ^G(k) y H(k) have 

j = o 
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Appendix E. Miscellanea on the Novel State Space Compensator 


E.l. An Example to Show the Compensation Principle of a State Space Filter 

Suppose there is a 2 nd -order linear time-invariant system, with a transfer function 
given by 


U (s) s 2 +a 1 s + a 0 


(E.l) 


f x = Ax + B« 

Its state equation and output equation are given by < , and the observable 

[ v = Cx + Du 


matrices are 


A = 


~ a \ 1 
-a 0 0 




0 

a n 


,C = [1 0],Z) = 0 


(E.2) 


Let x = [y y f , then the undelayed output is 


(E.3) 


With the state space compensation y = C<Px + C'FBw , the predicted output can be 
calculated by 

y p = 011*1 + 02*2 + ^12« 0 W (E.4) 

where <j> n and y/ r _ are elements of the matrix <h and its integration . Comparing Eqs. 
(E.3) and (E.4) gives y p =<t> n y+(/> l2 x 2 + y/ n a 0 u . It indicates that the continuous state space 
filter uses the system states and the control input to predict. The weighting factors of the 
relevant tenns are determined by the state transition matrix and its integral. Because the 
state variable y carries the velocity infonnation, the prediction makes use of the velocity. 
If the system is of 3 rd -order, the acceleration will be used to predict, and for higher order 
systems, more system information is required to make a prediction. 
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On the other hand, for a discrete system, say J 


n x z + n 0 
z 2 + m i z + m 0 


which is the 


discrete form of the system given by (E.l), the predicted output becomes 

y(k + I) = t/> l \y(k)- ( k - 1) + {n x y/ z u + n 0 y/~ u )u(k) + n 0 u (£-1) (E.5) 


Therefore, the discrete state space filter uses the current and past system information 
(states and input) to predict. 


E.2. The Filter Coefficients in Terms of the Eigenvalues 

A 2 nd -order system can be characterized by its natural frequency co n and damping 
ratio C ■ In terms of co„ and £ , the parameters in the transfer function in (E.l) can be 
written as 

o, = a 0 = <tf t (E.6) 
Let the two eigenvalues of this system be denoted by A, and A, , and the following 
relationships become straightforward 

a 1 =-(A l +A 2 ),a 2 =A l A 2 (E.7) 

As stated in Appendix A, the state transition matrix and its integral can be 
obtained in tenns of the system eigenvalues. The results are different for three cases with 
regard to the two eigenvalues, i.e., two different real eigenvalues, two repeated real 
eigenvalues, or a pair of complex conjugates. 

Two unequal real eigenvalues 

This case corresponds to an over-damped system, with g > 1 . The two eigenvalues 
are given by 

Kl ~ ~ m n £ ~ CO n \IC 2 “1 (E. 8) 
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Eq. (A. 9) becomes <P = a 0 +a { A . Without the details of derivation, a 0 and a, are 


given as follows 


a, =- 


A [ e / ^ d — 
e A ' td _ e KU 

A, -A 


A f d 


(E.9) 


where t d is the time delay. The state transition matrix and its integral are given, 
respectively, by 


O Atj 


K -K 


- A 1 e ktd -A, A? ( e K ' d - e ^ )' 
e Ku _ e ^i d A^e^ - Ai 1 e^ ,d 


a 0 - a, a, a, 


(E. 10) 


¥ 


_ e Au 
A, - A 2 


+(4 — 2X^)e 


^2 *d 


-1 


+ 1 


-A 2 e 


A<d 


A, -A 

A I +2A 2 e A, ' d 1 1 

A,-A 2 A 1 (Aj-A 2 ) A^(A l -A^) A, A 

l-or„ 


a, 


U x CI^OCq + Ct 2 OC\ 

6c n 1 


(E. 11) 


The D matrix of the compensated system can be calculated by 


( A 2 e*'' i +{Al 1 - 2/f ) e 

K 


D = C¥B = A l A 7 


= l-a n 


M l d \ 

— 1 


From (E.4), the prediction is given by 


(E.12) 


y P = 


Aifi^-A^tM* , A l A 2 (e^-e^' d ) 


A I -A 1 Aj—A^ 

= (a 0 -a l a)y + a l x 2 +(\-a Q )u 


A 2 e^ ld + (Aj -2A 2 )e 
A l —A 2 


A^td 


— 1 




(E-13) 
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For formulas (E.10)-(E.13), the last equations are in terms of a 0 and a,, the 
parameters of the transfer function, and a 0 and or, , which are given in (E.9). The 
remaining cases use the same convention. 

Two repeated real eigenvalues 

This case corresponds to a marginally damped system, with £ = 1 . The 
equationswhich correspond to Eqs. (E.8) through (E.13) for this case are given as follows 


A, 2 = A = -co n C 


f ^ Ttj , At., 

a 0 =e -t d e 


<D 


1 + At d t d 
-A 2 t d 1 -At d _ 


(E.14) 

(E.15) 


(E-16) 


T 


e d -At,e Md -1 


-e A,d + At d e A,J +1 

I 2 

, 2e A ‘ d — At d e A ‘ d —2 


A 


D = CTB = 1 - e A ' d (1 -At d ) 


y P =(l + At d )e A,d y + t d e A,i x 2 +[l-e A,d (\-At d )\i 


(E.17) 


(E.18) 

(E.19) 


Two complex conjugates 

This case corresponds to an under-damped system, with f <1 . The 
equationswhich correspond to Eqs. (E.8) through (E.13) for this case are given as follows 


A, 2 = ± i<°n J l -C 2 = a ±.iP 


(E.20) 


The relations given in (E.7) becomes 

a l =-2a,a 0 =a 2 +]3 2 (E.21) 
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a 0 = —^-{P sin pt d - a sin pt d ) 

gCCt d 

a\=—sinpt d 


<D 


P cos pt d + a sin pt d sin pt d 

- ( a 1 + p 2 ) sin pt d P cos pt d - a sin Pt d 


P 


r 


e a ' d sin Pt d 

P~ 


P + (asin pt d - p cos Pt d )e a 


(a 2 + P~)P 

(Pcos Pt d -asin Pt d )e a,d (/?' ~ a 1 } sin pt d + 2aP cos pt d ^ ^ 

P 1 [a 2 + P 2 )P e “ 20. 


(E.22) 


(E.23) 


(E.24) 


D = C*FB = 1 — 


a 


cos pt d - —sin Pt : 


(E.25) 




f a \ e a ‘ d 

cos Pt d -\ — sin pt d e°“ d y sin pt d x 2 + 

P ) P 


1- 


a 


cos pt d ——sin pt d \e a 


(E.26) 


For a 3 ld -order system 


IM. 

U(s) 


b 2 s 2 +b t s + b 0 

3 2 ? 

s a 2 s + ci^s + Uq 


its state transition matrix and its 


integral can be calculated by 


<S> = a 0 \ + a x A + a 2 A 2 


a 0 —a 2 a x +(#2 ~ a \ 

)a 2 

a , - a 2 a 2 

a 2 

—a x a x + (a,a, -a 0 ) 

1 a 2 

a Q - a, or. 

a, 

— a 0 or, + a x a 2 a 2 


-a 0 a 2 

«0 


or. 


1 -a n 


T = | « 0 — l—a x a 2 or, — (1 -a 0 ) + a 2 

a n 


-a 0 a 2 a 0 - 1 


a n 


(E.27) 


(E.28) 


Depending on the eigenvalues of the matrix A , the expressions of a 0 ,a x and a, may be 
different. Usually, a 3 rd -order aircraft model has one real and two complex conjugate 
eigenvalues. Denote them by A l =r,A 23 =m±n, then the following relations hold 
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(E.29) 


= — r(m 2 +n^,a 2 = 2 rm + (m 2 + n 2 a } = — (r + 2m) 


And the expressions of a 0 ,a t and a 2 are given by 




«0 - . 
A 


«o - . 
A 


-r(ra-a 2 + p 2 ) r(2a-r) , , , 

' e'” 1 ' sin fit d — e“'' ccw pt d + (A — r" + 2m)e” 


P 

r 2 - a 2 + p 2 


P 


P 


e a ‘ d sin /5t d - 2 ae“ d cos Pt d + 2ae" 


a-r 

~T 


e m ‘ sin pt d — e a ‘ d cos fit d + e" 


where A = r 2 + a 2 + p 2 - 2 ra . The matrix D of the compensated system is 


D = CTB = b 2 (a 1 - a 2 a 2 ) + b x a 2 + b, 


l-a n 


(E.30) 


(E.31) 
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Appendix F. Calculation of the Accelerations in the Topodetic Frame 

The topodetic frame is also called the north-east-up (NEU) frame. This is the 
frame in which most visual images in a flight simulator are generated. The six degrees of 
freedom (DOF) are: roll, pitch, yaw (i.e., the three Euler angles), altitude, longitude and 
latitude. The displacements, velocities and accelerations in the 6-DOF of the topodetic 
frame are usually calculated from those in the 6-DOF of the body frame, which are 
available from the flight dynamics and equations of motion (EOM). Because the 
accelerations in the topodetic frame are not used in the visual image generation, they are 
not usually available from the simulation computer. As mentioned in Chapter 4 (4.5.1), 
they are necessary for applying the state space predictor in a flight simulator to 
compensate the transport delay in the visual system. The formulas to calculate the 
accelerations in the 6-DOF of the topodetic frame are derived as follows. These six 
accelerations can be divided into two groups — three angular accelerations (2 nd derivatives 
of the Euler angles) and three translational accelerations. 

F.l. Angular Accelerations 

The transformation of the angular velocities from the body frame to the earth 
frame is given by the kinematic equations 

<j) 1 tan 0 sin (j) tan 0 cos (j) p 

6 = 0 cos(j) - sin (j) q . (F.l) 

!// 0 sin (j) / cos 0 cos (j) / cos 0 r 

Denote v E K =\J> 6 f/] (velocity vector in the earth or topodetic frame for the 

kinematic equations), \ B K =\p q rf (velocity vector in the body frame for the 
kinematic equations), and 
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where the nine elements are given in Table F.l. 

Investigation shows that the second term of Eq. (F.4) is much smaller than the 
first term (the former is about one thousandth of the latter), and hence Eq. (F.4) may be 
approximated by discarding the second term, that is 

(F.7) 
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or 


0 


“l 

tan 6 sin (j) 

tan 0 cos (f) 

P 

0 

~ 

0 

COS(j) 

- sin (j) 

q 

V 


0 

sin (j) / cos 0 

cos (j) / cos 9 

r 


This is exactly Eq. (4.34). 

Table F.l. Expressions of the nine elements of matrix Pg 2£ Q 


Element 

Expression 

Ctji 

v (cos <j)sin 9 cos y/+sin <j)sin i//) + w(co5 tpsin i//— sin (p sin 6 cos i//) 

a n 

— u sin 9 cos y/ + v sin (p cos 9 cos y/ + wcostpcos 9 cos y/ 

OC u 

—u cos 9 sin y/ — v ( sin tpsin 9 sin y/ + cos tpcos y/) + w[sin tpcos y/ - cos <p sin 9 sin y/) 

«2I 

v (cos <p sin 9 sin yr - sin tpcos y/)— w(cos tp cos y/ + sin tpsin 9 sin y / ) 

^22 

—u sin 9 sin y/ + v sin tpcos 9 sin y/ + wcostpcos 9 sin y/ 

oc 23 

u cos 9 cos y/+v ( sin tpsin 9 cos y/ - cos tp sin y/) + w^sin tp sin y/ + cos tp sin 9 cos y / ) 

®31 

-v cos tp cos 9 + wsintp cos 9 

OC 32 

u cos 9 + v sin tpsin 9 + w cos tpsin 9 

a n 

0 


F.2. Translational Accelerations 


The transformation of the angular velocities from the body frame to the earth frame is 
given by the navigation equations 


T 


i 

h 

— 


cos 9 cos iff sin (/) sin 9 cos iff - cos (j) sin (// cos (f) sin 9 cos y/ + sin tpsin (// 

cos 9 sin (// sin (p sin 9 sin (// + cos (p cos (// cos (p sin 9 sin (// - sin (p cos (// 

sin 9 - sin (p cos 9 - cos (p cos 9 


(F.9) 


Denote \ E N = [A / h (velocity vector in the earth or topodetic frame for the 


navigation equations) and \ b n = \ii v w T (velocity vector in the body frame for the 
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navigation equations). Denote the transformation matrix for the translational velocities 
from the body frame to the topodetic frame as 


P r = 

1 B2E 


cos 0 cos y/ 
cos 6 sin (// 
sin 0 


sin (psin 0 cos (// - cos (j) sin (// 
sin (p sin 6 sin y/ + cos (p cos y/ 
- sin (p cos 6 


cos (psin 0 cos y/ + sin (p sin y/ 
cos (psin 9 sin (// - sin (p cos (// , 
- cos (p cos 6 


(F. 10) 


then Eq. (F.9) can be rewritten in a compact fonn 

\ E = p r y B 
y N l B2E y N 4 


(F -11) 


The translational accelerations may be calculated by taking the derivatives on both sides 


of Eq. (F.l 1), and the final result is 


11 


u 



/ 

= p r 

A B2E 

V 

+ T 

0 

h 


w 


¥ 


(F.12) 


This is the same as Eq. (4.35), and the matrix T is available in Table 4.2. 

In the VMS, the longitude and latitude are expressed in angles, therefore, the 
tenns X and l have to be divided by the corresponding radii to change to angles. The 
results are 


h + 

a 2 b 2 

V 

| ~ b 2 +(ff 2 -b 2 ^cos 2 /lj 


1 cos A 


a , 

+ /? 


cos 


4 + .< 




V a 

with a=2.092565e+7 ft, the earth equatorial radius and b = a 
where f=l/298.257 the earth flattening parameter. 


(F. 13) 


(F.14) 


(1-/) , the earth polar radius, 
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